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We present a theory of the hnear response of a Bose-condensed gas to external perturbations at 
finite temperature. The theory developed here is the basis of a recent quantitative explanation of the 
measurements of condensate excitations and decay rates made at JILA [D. S. Jin et al, Phys. Rev. 
Lett. 78, 764 (1997)]. The formalism is based on a dynamic, number-conserving, mean-field scheme 
and is valid in the collisionless limit of well-defined quasiparticles. The theory is gapless, consistent 
with the generalized Kohn theorem for the dipole modes, and includes the time-dependent normal 
and anomalous averages, Beliaev and Landau processes, and all relevant finite size effects. The 
important physical process where the thermal cloud is driven directly by the external perturbation 
is explicitly included. This is required for consistency with the dipole modes and is also needed to 
explain the JILA observations. 
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I. INTRODUCTION 

Some of the earliest experiments on Bose-Einstein con- 
densates (BEC) in trapped gases involved measuring 
the energies and decay rates of condensate excitations 
P, 0, IM^' These measurements provide a sensitive 
probe of condensate dynamics and a unique opportu- 
nity for quantitative tests of thermal field theories. For 
this reason condensate excitations continue to attract 
widespread interest, and a variety of new measurements 
have been made, notably for scissors mode oscillations 
la , in highly elongated traps (^1, via Bragg spectroscopy 
0, and in vortex condensates [g]- In this paper we de- 
rive a second order quantum field theory which can be 
applied to these experiments. This theory has recently 
been used to obtain good agreement with finite tempera- 
ture measurements made at JILA in 1997 0, which have 
been the subject of much discussion in the field 

Theoretical calculations of condensate excitations in 
dilute Bose gases are often based on the quasiparticle 
description introduced by Bogoliubov |1C| . However, ac- 
curate calculations beyond the Bogoliubov theory are dif- 
ficult for a variety of reasons. First, the various higher 
order effects which occur have to be treated consistently 
to avoid problems of infra-red divergences and to obtain 
a gapless spectrum, as required by the Hugenholtz-Pines 
theorem Second, one must in principle deal simulta- 
neously with the strong modification of low-energy states 
due to interactions and the substantial single-particle ef- 
fects that are present at finite temperature. Third, a dy- 
namic description of both condensed and non-condensed 
atoms and their mutual interaction is required. While 
the condensate dynamics is well described by a sin- 
gle nonlinear equation (the Gross-Pitaevskii equation or 
GPE 13) J evolving the non-condensate is a much more 
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complicated numerical problem. Finally, experiments in 
trapped gases may have substantial finite size effects so 
that approximations valid in the thermodynamic limit 
may not be appropriate. 

The theory we develop in this paper addresses all these 
issues using a full quasiparticle description of the non- 
condensate. The dynamic coupling between condensed 
and uncondensed atoms is described using perturbation 
theory. It has been stated in the literature that a per- 
turbative calculation can not explain the JILA measure- 
ments 0, Q| , but our recent calculations show that this 
is not the case 0, llSj . A perturbative approach is per- 
fectly adequate, although it must be applied with great 
care. It is particularly important to include the effect of 
the external perturbation on the non-condensate, as well 
as condensate phase and number fluctuations and vari- 
ous finite size effects, all of which are relevant in the JILA 
experiment. An advantage of the perturbative approach 
is that explicit expressions can be obtained for various 
response functions and this simplifies numerical calcu- 
lations, allowing the effects of low-energy phonons and 
high-energy single-particles to be included in the same 
calculation. 

The study of the dilute Bose gas has a lon g h istory 
dating back to the work of Bogoliubov in 1947 nul, who 
described the excitations in terms of weakly-interacting 
quasiparticles. Beliaev extended the Bogoliubov ap- 
proach by including interactions between quasiparticles 
using perturbation theory [l^ . This approach was fur- 
ther developed by Hugenholtz and Pines 11] and by 
Popov and Fadeev [TtI fl^ who applied it to calculations 
at finite temperature. The Bogoliubov method was ap- 
plied to inhomogeneous systems by de Gennes and Fetter 
19.. .20] , while the Beliaev-Popov formalism has recently 
been reexamined by Shi and Griffin 21] and extended 
to trapped gases by Fedichev and Shlyapnikov The 
excitation spectrum at zero and finite temperature was 
also studied by Mohling, Sirlin and Morita using many- 
body perturbation theory These calculations have 
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recently been extended to trapped gases by Morgan . 

Dynamical properties of homogeneous and trapped 
condensates are well-described at low temperatures by 
the GPE Condensate excitations can be stud- 

ied by applying linear response theory to the GPE 
I25I l26l| . reproducing the Bogohubov quasiparticle re- 
sults. Measurements of excitations at low-temperature 
are in good agreement with predictions based on this 
approach P, l2a. l27j . The GPE can be generalized to 
include the dynamic coupling to non-condensed atoms, 
which is important at finite temperature. Linear re- 
sponse theoryhas been applied to the generalized GPE 
by Giorgini ^§\, and leads to results which are equiva- 
lent to the Bcliacv-Popov formalism applied to trapped 
gases. Linear response theory at finite temperature has 
also been discussed by Minguzzi and Tost [23, Guilleu- 
mas and Pitaevskii jsfl and by Rusch and Burnett [sijl . 
while the related dielectric formalism has been applied in 
this context by Bene and Szepfalusy and Reidl et al. 
[T3II . Recently, the GPE has been used to study the finite 
temperature dynamics of highly-occupied modes in Bose 
gases JJSj jMj J and a stochastic GPE has been developed 
to include the effect of quantum fluctuations |3a, [Sy, |33 ■ 

Although low-temperature excitations can be stud- 
ied straightforwardly using the GPE, finite-temperature 
measurements, particularly those made at JILA 0], have 
proved much harder to explain. The starting point for 
many recent calculations of condensate excitations at fi- 
nite temperature is the Hartree-Fock-Bogoliubov (HFB) 
theory, which has been summarized by Grifhn (38i |. The 
full form of this theory is not gapless, however, and the 
HFB theory is therefore generally used in an approxima- 
tion, often called the Popov approximation (HFBP), in 
which the anomalous (pair) average is neglected. The 
HFBP theory has been ap plied to trapped gases by a 
number of authors |3^E3k41J, and a gapless extension 
of the HFB formalism which includes the anomalous av- 
erage has also been developed (GHFB) [12, The 
various HFB approaches have been summarized recently 
in 

The application of the HFB approach to the 1997 JILA 
experiment 2] is particularly interesting. In this exper- 
iment the energies of the lowest-energy modes with ax- 
ial angular momentum quantum numbers m = 2 and 
m — were measured as a function of reduced temper- 
ature t = T/T^, where T is the absolute temperature 
and is the BEG critical temperature for an ideal gas. 
The m — 2 mode was observed to shift downwards with 
t, while the m = mode underwent a sharp increase in 
energy at t 0.6 towards the result expected in the non- 
interacting limit. HFBP calculations give good agree- 
ment with the JILA data for low temperatures but can 
not explain the results for t > 0.6 ^^]. Good agreement 
for all t for the 171 — 2 mode was obtained using the 
GHFB theory, and also within the dielectric formalism 
[l3l |. However, neither method was able to explain the 
upward shift of the m = mode, and an analytical cal- 
culation based on linearizing the generalized GPE also 



predicted downward shifts for both modes '28|. 

A possible explanation for the behaviour of the m = 
mode has been given by Bijlsma, Al Khawaja and Stoof 
[is !. . They emphasized the importance of the relative 
phase of condensate and non-condensate fluctuations and 
showed that the experimental results for m ~ can be 
qualitatively explained by a shift from out-of-phase to in- 
phase oscillations at high temperature. Recently, Jack- 
son and Zaremba obtained good agreement with the 
JILA experiment for both modes and also for the Ox- 
ford scissors mode experiment Their calculations 
are based on the Zaremba-Nikuni-Griffin formalism [4^ 
and use a generalized GPE for the condensate coupled to 
a non-condensate modelled by a semiclassical Boltzmann 
equation. Despite its successes, however, this approach 
neglects the phonon character of low-energy excitations 
as well as the anomalous average and all Beliaev pro- 
cesses. It is not a priori obvious that these are valid ap- 
proximations. For example, in a weakly-interacting Bose 
condensate the anomalous average is as large as the con- 
tribution to the normal average from interactions |24| . 
and it can be very significant near a Feshbach resonance 
[4^ Isof . The good agreement with the JILA results for 
the m = 2 mode within the GHFB theory Q also sug- 
gests that the anomalous average can be important. Fur- 
thermore, Beliaev processes have been dir ectly observed 
in a number of recent experiments (Sj, l52l Issj . It is 
therefore important to develop a theory which includes 
these effects in a consistent manner, and we describe such 
a theory in this paper. An alternative kinetic theory 
which also includes these processes consistently has been 
developed by Walser et al. and Wachter et al. 

As well as measuring the energies of excitations as 
a function of temperature, experiments can also study 
the decay rates. The HFB theories do not describe the 
damping of excitations, however, because they treat the 
thermal cloud statically, whereas damping arises from 
the dynamical interaction between condensed and un- 
condensed atoms which removes energy from the con- 
densate fluctuation. Damping of excitations has been 
studied theoretically by a number of authors (see for ex- 
ample [IlllllllpflillilililiiiS among many 
others) and there is generally reasonable agreement with 
experimental measurements. 



A. Outline of the paper 

In this paper we present a theoretical description of 
condensate excitations which has recently been used suc- 
cessfully to explain the JILA experimental results for 
both the m = 2 and m — Q modes 0. Following 
Giorgini 2^ , we use linear response theory with a gener- 
alized GPE. We obtain this equation using the number- 
conserving formalism of Gastin and Dum j6lj with the 
result that it differs in detail from the more usual form 
obtained via spontaneous symmetry breaking. The lin- 
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ear response treatment has the advantage that it closely 
models the experimental procedure where excitations are 
created by small modulations of the trap frequencies. In 
particular we can easily include condensate phase and 
number fluctuations and the effect of the external per- 
turbation on the non-condensate, which were neglected 
in our earlier time-independent approach |24ll62j |. 

This paper is organised as follows. In Sec. |n] we 
summarize the Castin-Dum number-conserving formal- 
ism and obtain equations of motion for the condensate 
and non-condensate. We also introduce the quasiparti- 
cle basis and discuss its relevant properties. In Sec. IIIII 
we briefly discuss the UV renormalization of the the- 
ory and outline its justification. In Sec. IIVI we develop 
the linear response formalism for a Bose gas, starting 
in Sec. IIV Al with the simplest approach based on the 
GPE and then considering in Sec. IIV Bl the higher order 
corrections which arise from the generalized GPE. Vari- 
ous issues related to the resulting theory are discussed in 
Sec. 13 In particular, in Sec. IV Al we consider the homo- 
geneous limit which can be used to show that the theory 
is gapless and to obtain the conditions for its validity. 
In Sec. IV Bl we show that non-condensate collisions are 
of higher order in the small parameter of the theory and 
can therefore be neglected for the purposes of this work. 
In Sec. IV Cl we show that the dipole modes are obtained 
correctly to within the order of the calculation. The ap- 
pendices contain further details of the calculations. 

We briefly comment on our notation, which is poten- 
tially complicated as a number of similar yet distinct 
quantities are required in the theory. We denote by a 
subscript '0' or '2' the relevant approximation to the con- 
densate wavefunction, '0' for the GPE and '2' for the 
generalized GPE. Any time or frequency dependence is 
written explicitly (unless stated otherwise), so if none is 
given the corresponding quantity is calculated in equi- 
librium. Spatial dependence is usually written explicitly 
but is occasionally neglected when it is obvious. Most 
quantities in the theory depend on temperature, but to 
avoid overloading the notation we do not indicate this 
explicitly. 



II. SUMMARY OF TIME-DEPENDENT 
NUMBER-CONSERVING FORMALISM 

Theoretical descriptions of dilute Bose condensates 
usually start from the assumption, first applied in this 
context by Bogoliubov ^Td\, that the U{1) gauge symme- 
try of the Hamiltonian is spontaneously broken. While 
this procedure is undoubtedly useful there are some tech- 
nical concerns with the violation of particle number- 
conservation that is involved. In fact, it has been shown 
by Girardeau and Arnowitt (g^I, Gardiner |63|, and 
Castin and Dum '6l'| that symmetry breaking is not nec- 
essary and the Bogoliubov theory and higher order cal- 
culations can be obtained via a number-conserving ap- 
proach. In this paper we use the time-dependent for- 



malism of Castin and Dum [6l||, with minor extensions 
to finite temperature, to obtain the equations of motion 
for the condensate and non-condensate. At leading or- 
der there is no difference between this approach and a 
symmetry-breaking formalism. At higher order, however, 
the number-conserving approach describes various finite 
size effects which are not contained in symmetry-breaking 
calculations. We found in J3J that these effects are impor- 
tant for explaining the JILA experimental results within 
this formalism, although it is an open question to what 
extent they are significant more generally. 

We have chosen to use an explicitly number-conserving 
approach rather than a symmetry-breaking treatment be- 
cause we will make regular use of the Bogoliubov quasi- 
particle basis to expand various wavefunctions. It is 
therefore important to have a well-defined, complete ba- 
sis set and in particular to deal correctly with the zero- 
energy condensate modes which appear. This is simplest 
in the number-conserving approach where these states 
occur naturally and are treated on the same footing as 
those with non-zero energy. In particular, there is no 
'missing eigenvector' in this approach and issues related 
to the orthogonality of excitations to the condensate are 
dealt with explicitly. 

In this section we summarize the Castin-Dum number- 
conserving formalism. The original derivation given in 
|6l| was restricted to near zero temperatures because it 
was assumed that the condensate population Nq was ap- 
proximately equal to the fixed total number of atoms N. 
The changes required to adapt the approach to pertur- 
bative calculations at finite temperature are described in 
more detail elsewhere [6^. In this section we quote the 
equations that result and discuss the properties that are 
relevant to the current work. 



A. Equations of motion 

We start from the usual Heisenberg equation of motion 
for the field operator ^'(r, t) which annihilates a particle 
at point r and time t 

z;i^*(r,t) =i?,p(r,t)*(r,i) 

+ C/o*^(r,t)^'(r,t)*(r,t), (1) 

where i/sp(r,<) = -{h'^/2m)V^ + V{r,t) is the single- 
particle Hamiltonian containing the kinetic energy and 
any external potential. For simplicity we have as- 
sumed binary interactions between the particles char- 
acterized by an energy-independent contact potential, 
14in(r — r') = Uo6{r — r'), where Uq = Anfi^as/m and 
as is the s-wave scattering length. This is the standard 
approximation for three-dimensional, dilute Bose gases. 
As is well known, however, it leads to ultra-violet diver- 
gences which must be removed by renormalizing various 
quantities which appear in the subsequent development 
of the theory. This procedure is well understood and has 
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been discussed by a number of authors (see for example 
113, Ell among many others). It can be rigor- 
ously justified and we give a brief outline of the relevant 
arguments in Sec. IIIII 

Following the definition of Penrose and Onsager , a 
condensate is said to be present if the one-body density 
operator pi{r,r\t) — {^''^ {r' ,t)'^{r,t)) has an eigenvec- 
tor with a macroscopic eigenvalue much larger than all 
others. Here (. . .) denotes a statistical average in the 
appropriate ensemble. The exact condensate wavefunc- 
tion <i>ox(r, t) (normalized to one) and population No{t) 
therefore satisfy the equation 

y"dr'pi(r,r',t)$ex(r',t) = A^o(i)*cx(r,t), (2) 

with NQ{t) ^ 1. We can write the full field operator in 
terms of condensate and non-condensate operators as 



*(r,i) -$ex(r,i)ao(i)+'5*(r,i), 



(3) 



where the condensate annihilation operator ao{t) is de- 
fined by projecting ^{r,t) onto $ox(r, i), while (55'(r,t) 
is obtained via the orthogonal projection 



ao{t) - ydr$:,(r,t)*(r,t), (4) 
J#(r,i)= /dr'Qex(r,r',i)*(r',0 = Qcx(i)*(r,t). (5) 



The orthogonal projector Qox(r,r',i) is defined by 

Qex(r, r', t) = Sir - r') - <i>ex(r, t)<f ^(r', i), (6) 

and the final part of Eq. Q defines a convenient short- 
hand notation for the action of this operator. The mean 
condensate population is given as usual by iVo(i) = 
{al{t)ao{t)). 

Eq. ^ leads directly to 



(aS(t)5*(r,t)) =0, 



(7) 



which is the defining relation from which equations for 
<l'cx(r, t) and 5^(r, t) follow 61]. The equation of motion 
for aj^{t)5^{r,t) can be written as 



al[t)5i>{v,t) 



dr'J2Mr,r',t), (8) 



k=0 



where the quantities Rk{r,r',t) are operators containing 
k factors of the non-condensate field 6^. In the presence 
of a condensate we can therefore assume that the effect of 
Rk+i is small compared to Rk and use this as a means of 
developing a systematic approximation scheme. Includ- 
ing progressively more of the Rk leads to a sequence of 
improving approximations to the exact condensate wave- 
function, which we denote by <i>o(r, i), <i>i(r, i), <i>2(r,t) 
etc. 



The leading order theory is obtained by neglecting all 
terms except Rq, which gives an approximate equation 
of motion for the condensate wavefunction ^'ex(r, t) 
$0(1", t) which is just the usual time-dependent GPE 



in^$o(r,t) 



+ NoUo\Mr,t)\Ho{r,t). (9) 



Here $o(r,t) is normalized to one and Xo{t) is an arbi- 
trary, real function of time which can be chosen at will 
to adjust the global phase of the condensate wavefunc- 
tion. At this level of approximation the mean condensate 
population Nq is independent of time. 

Improved approximations for the condensate wave- 
function can be obtained by including more of the opera- 
tors Rk. It turns out that (-Ri) = 0, so <i>i(r, t) = $o(r, t). 
The next approximation therefore comes from including 
the effect of R2, $cx(i", t) 'i'2(r, t). The resulting equa- 
tion of motion is the generalized GPE (GGPE) 



ift^$2(r,i) 



i?sp(r,t)-A2(t)J $2(r,i) 

+ {No{t) + ANo)Uo\^2ir,t)\H2{r,t) 
+ 2C/on(r, i)$2(r, t) + C/om(r, t)$*(r, t) 
-f{r,t), (10) 



where the function ,f{r,t) is defined by 



/(r,i) = y"dr'[/o|$2(r',i)p[$2(r',t)n(r,r',t) 

+ $^(r',t)m(r,r',i) , (11) 



and AiVo is defined by 



ANo 



(No)' 



No 



1. 



(12) 



The contribution to the GGPE from /(r,t) was first 
obtained by Castin and Dum 61] and does not appear 
in symmetry-breaking approaches to the theory of BEG. 
We will see later that it is related to the strict enforce- 
ment of orthogonality between the condensate and non- 
condensate. The contribution from statistical fluctua- 
tions in the condensate population (AiVo) is new and 
comes from the extension of the Gastin-Dum analysis to 
finite temperature . In practice, these fluctuations are 
negligible, both because they are intrinsically small even 
for condensates containing only a few thousand atoms, 
and also because the properties of the condensate are in- 
sensitive to the exact number of atoms when Nq is large. 
For completeness, however, we retain this term in the 
GGPE during the development of the theory. An ex- 
pression for ATVo in terms of Bogoliubov quasiparticles 
is given in Appendix IdI 

The off-diagonal non-condensate density matrix 
n{r, r', t) and anomalous average m(r, r', t) which appear 
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in Eqs. ifTUI) and lfTT|l are defined by 

n(r,r',i)-(At(r',t)A(r,t)), (13) 
m(r,r',i) = (A(r',t)A(r,t)\, (14) 



where A(r, t) is the number-conserving, non-condensate 
field operator IHl El m 



A{r,t) = ^=al{t)5i>{v,t). 



(15) 



For the diagonal parts of the non-condensate density and 
anomalous average we use the notation n{r, t) = n(r, r, t) 
and TO(r, t) = m(r, r, t). 

At this level of approximation, the mean conden- 
sate population No{t) is no longer independent of time. 
The condensate population is determined by the non- 
condensate population Nnc{t) = Jdrh{r,t) because the 
total atom number is fixed 



No{t) + N,,,{t) = N - 



dNo 
dt 



dt 



(16) 



In fact we use this relation to determine the condensate 
population in terms of N and the temperature T. The 
parameter Nglt) is therefore just a shorthand notation for 
N — Nnc{t)- Equation l|16|l shows that the evolution of 
Nq (t) depends quadratically on the non-condensate field 
operator. This means that its time-dependence is only 
relevant when we consider the effect of i?2 on the con- 
densate wavefunction, i.e. in the GGPE of Eq. (|10|l . and 
at leading order we can assume that iVo is independent 
of time. 

As a consequence of the time-dependence of No{t), the 
parameter A2 {t) in Eq. H1U|) has an imaginary part (t) 
in order to keep the normalization of the wavefunction 
constant, /dr |$2(r,i)P = 1- This condition gives 



A^2W = ^ Jdr [Uorhir,t)^fir,t) 



c.c. 



h dNo 

2Nn dt ' 



(17) 



where c.c. stands for 'complex conjugate' and the last 
line follows from the equation of motion for NQ{t) to the 
order of this calculation. The imaginary part of X2{t) 
is therefore determined by the equations of motion and 
only the real part is arbitrary (as before, this corresponds 
to choosing the global phase of the condensate wavefunc- 
tion). Eq. H17|) shows that at this level of approximation, 
the dominant processes changing the condensate popula- 
tion involve pair excitation into and out of the conden- 
sate. 

The GGPE must be combined with an equation of mo- 
tion for the non-condensate field operator A{r,t). Cor- 
rect to the order of this calculation, this is given by 



7,1 
dt VAt(r,i) 



A(r,t) 



= C2{r,t) 



A{r,t) 
At(r,t) 



(18) 



where the matrix operator £2(1", t) is defined by 



C2ir,t) 



L2(r,t) M2{r,t) 
-m{v,t) -L*^{v,t) 



(19) 



i2(r,i) = -ffsp(r,t) - A2(t) +iVof/o|$2(r,t)|' 

+A^oC/oQ2(t)|*2(r,t)pQ2(i), (20) 

M2(r, t) = NaUoQ2ml{^, t)Q*2(t)- (21) 
The notation Q2{t)\^2{'^,t)\'^Q2{t) means, for example. 



Q2(i)|$2(r,t)pQ2(t) 



A(r,i) 



JJd^r'dh" Q2(r,r',i)|<f2(r',i)|2 

xQ2(r',r",i)A(r",t), (22) 

where (32(1", r', t) is defined as in Eq. ® with $cx(r, t) 
$2(1", i). Equation 118|) differs from the form one would 
obtain in a symmetry-breaking approach by the appear- 
ance of the orthogonal projectors Q2- These preserve or- 
thogonality between the time-dependent condensate and 
non-condensate and are unique to the number-conserving 
formalism. A subtle issue concerns the appearance of 
$2(r, i) rather than $o(r,i) in Eqs. This is 

necessary to obtain sensible expressions for energy shifts 
and is discussed further in Appendix IbI 

The GPE of Eq. ®, its generalization of Eq. ^ 
and Eq. IjlSI) for the non-condensate (together with the 
definitions of all relevant quantities) define the time- 
dependent, number-conserving, mean- field theory we use 
in this paper. Before applying these equations, however, 
we first describe the quasiparticle basis obtained from the 
eigenstate solutions of Eq. H18|l and discuss some of its 
properties. 



B. Bogoliubov Quasiparticles 

Castin and Dum 61] have shown that the solutions of 
Eq. H18|l can be written in terms of a basis expansion 
with time-independent quasiparticle annihilation opera- 
tors and time-dependent eigenmodes. This is very use- 
ful in the subsequent development of the theory. The 
static case defines a convenient basis for the calcula- 
tion while the dynamic case gives expressions for the full 
off-diagonal time-dependent non-condensate density and 
anomalous average. 



1. Static case 

We consider first the situation that the condensate 
wavefunction satisfies the static limit of the ordinary 
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GPE of Eq. . If the external potential is independent 
of time, we can find solutions $o(r,i) = $o(r) which 
satisfy 



Hspir) - Ao + iVoC/o|$o(r)P $o(r) = 



(23) 



where Aq (the condensate eigenvalue) is a constant, ap- 
proximately equal to the chemical potential of the sys- 
tem, and the single-particle Hamiltonian iJsp(r) contains 
the kinetic energy and any static trap potential 



2m 



+ ^4rap(r). 



(24) 



We now consider the matrix operator £2(1", t) for this 
static case with the replacements $2(r,t) '^'o(r), 
A2(i) = Aq. This gives the operator £o(r), defined by 
making these changes in Eqs. H19|l - H21|) 



^°^^)-U>o(r) -Ll{v) 



(25) 



is an eigenvector with eigenvalue — and norm —Xi 
(cti = is the first Pauli matrix). If $o(r) is a sta- 

ble ground state solution of the GPE, then states with 
positive norm have positive or zero energy while states 
with negative norm have negative or zero energy. This is 
not true, however, if <I'o(r) contains a vortex. In this case 
there exist anomalous modes with positive norm and neg- 
ative energy which have been the source of much recent 
discussion (see for example 71, 72, 73, 74, 75j). 

The operator £0 (r) is diagonalizable so its eigenvectors 
form a complete set and there are no missing modes in 
this formalism. In particular, there are two zero-energy 
solutions to Eq. H28|l which we denote by ±0, according as 
their norm is ±1. These two states involve the condensate 
wavefunction and are given by 



A'+o(r) 
A'-o(r) 



w+o(r) 
V-o{r) 



$o(r) 






(32) 
(33) 



Lo(r) = i?sp(r) - Ao + N^U^\<^^{y)\' 

+Qo|$o(r)|'Qo, (26) 

Mo(r) = iVoC/oQo*o(r)Qo- (27) 

The eigenstates and eigenvalues of £0 (r) are the solutions 
of the Bogoliubov-de Gennes (BdG) equations 



(28) 



where the two-component spinor Xi (r) is defined in terms 
of quasiparticle wavefunctions Uiiv) and Vi{r) by 



(29) 



and the are the quasiparticle energies. These solutions 
define the static quasiparticle basis we use in this pa- 
per. At this level of approximation, the theory is simply 
the Bogoliubov theory with orthogonalized quasiparticle 
wavef unct ions . 

The eigenvectors Xi{v) have a number of important 
properties that we will need later |^ [t^- They obey 
a generalized orthonormality relation which can be sum- 
marized by 



jdvX}{v)a:,X,{v)^ 



(30) 



1 



where Xi = il is the norm of Xi(r) and (T3 

is the third Pauli matrix. If Xi(r) is an eigenvector with 
eigenvalue and norm Xi then 



A'_,(r) = aiA;*(r) 



(31) 



An advantage of the number-conserving formalism, is 
that these two condensate modes obey the same rela- 
tions as states with non-zero energy and can therefore be 
treated on exactly the same footing. In particular, they 
are related by the same interchange of u and v functions 
that applies to other pairs of positive/negative energy 
states and they obey the orthonormality condition, which 
in this case simply expresses the orthogonality of modes 
with non-zero energy to the condensate 

Jdrmr)u,ir)^ Jdr <Po{r)v,{r) ^ 0, {i^±0). (34) 

It turns out to be convenient to form a linear com- 
bination of the two condensate modes to obtain vectors 
describing changes in the condensate norm and phase 

These states satisfy the generalized orthogonality of 
Eq. H3U|) with all other states but have zero norm. They 
can be projected out of any expansion using the result 



i ldrXl^{r)<jMr) = l. 



(36) 



The completeness of the basis vectors Xi{r) gives the 
decomposition of unity 



1 
1 



(37) 



where the sum is over all states. The fact that these vec- 
tors form a complete set of states is very important for 
our subsequent use of them as a basis to expand conden- 
sate and non-condensate fiuctuations. 
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2. Dynamic case 

In the dynamic case, the non-condensate field opera- 
tor A(r, t) can be written in terms of time-independent 
quasiparticle creation and annihilation operators (respec- 
tively /3|, and time-dependent eigenmodes 1(5 1| . Since 
A(r,t) is orthogonal to the condensate we can write 



(38) 



i>0 



where the summation is over all non-condensate modes 
with positive norm, Xi = +1- The spinors <%i(r, t) are de- 
fined in terms of time-dependent Bogoliubov wavefunc- 
tions Ui{r,t) and Vi{r,t) by 



Ui{r, t) 
Vi{r, t) 



u*{r,t) 



and evolve according to 
d 



ih—X,{v,t)=C2{r,t)X,{r,t). 



(39) 



(40) 



In the equilibrium case, the spinors Ai(r,t) are related 
to the static basis functions defined above by Xi (r, t) = 
Xi{r:)e-"-^^/^. 

In this formulation, the quasiparticle creation and an- 
nihilation operators are independent of time 



dt 



dt 



= 0. 



(41) 



It follows from the orthonormality relations and the def- 
inition of A(r,t) that the quasiparticle operators obey 
bosonic commutation relations (neglecting the subspace 
with TVq = 0). We note that Eq. gUJ) [with $2(r,i) 
$o(r, i) in £2(1", t)] also gives the evolution of the con- 
densate to leading order because in this case it reduces 
to the GPE of Eq. © . This ensures preservation of or- 
thogonality between the condensate and excited states to 
the order of the calculation in the time-dependent case. 

Equations and (|in|l provide the full equation of 
motion for the non-condensate from which n{r, r', t) and 
m{r,r',t) can be constructed. A representation of these 
quantities in terms of the spinors Xi{r,t) is given in Ap- 
pendix |X3| Eqs. ljll7l) - (|ll9|) and is very convenient for 
calculations. They can also be written in a compact 
form using the time-dependent quasiparticle wavefunc- 
tions Ui{r,t) and Vi{r,t) as 



i(r,r',t) = ^u,(r,t)<(r',i)A^, 



i>0 



+v:{r,t)v,{r',t){N, + l), (42) 



i(r,r',t) = ^u,(r,t)<(r',t)iV, 



i>0 



+<(r,t)M,(r',t)(iV, + l), (43) 



where we have assumed that (PjPj) = NiSij and {PiPj) = 
0, which are appropriate for non-interacting quasiparti- 
cles. As a consequence of Eq. H41|l . the quasiparticle pop- 
ulations {Ni\ are independent of time, as has also been 
shown by Walser et al. j^. Substituting Eqs. and 
(|^ into Eq. ((111) gives the form of /(r, t) quoted in ^1 



/(^'^) ^ ^y,c:{t)N,u,{r,t) + cS){N, + l)<(r,i), 

(44) 

where the scalars Ciit) are defined by 

c,{t) = jdvN^Uo\^2(r,t)\''[^*{v,t)u,{Y,t) 

+ $2(r,t)w,(r,i)]. (45) 

We will be interested in the case that the quasiparticle 
populations {Ni\ correspond to thermal equilibrium so 
that they are given by the Bose-Einstcin distribution 



1 



g/3(ei-[p-Ao]) _ 1' 



(46) 



where is the chemical potential, (3 = and fee is 

Boltzmann's constant. The difference between ^ and Aq 
is of order 1/Nq, but for systems with a finite number of 
particles this difference must be accounted for otherwise 
it is impossible to find equilibrium solutions that satisfy 
the constraint on N near the critical temperature (and 
beyond). We note that most quantities in the theory de- 
pend on temperature via their dependence on these pop- 
ulations [either explicitly as in n{r, r', t) and m(r, r', t) or 
implicitly as in <I>2(r, i) and NQ{t)]. 



III. ULTRA-VIOLET RENORMALIZATION 

The GGPE of Eq. (fTUIl contains the diagonal part of 
the anomalous average m{r, t) which is ultra-violet (UV) 
divergent. The reason and cure for this problem are well- 
known and we give a brief summary of the argument in 
this section. 

The UV divergence in the anomalous average arises be- 
cause of the use of the contact potential approximation 
to describe particle interactions. A truly ab initio theory 
would have to start by describing these interactions using 
the actual interatomic potential. In contrast, the contact 
potential is the zero-energy, zero-momentum limit of the 
homogeneous, two-body T-matrix describing the scatter- 
ing of two particles in a vacuum. Detailed arguments 
which show formally how the bare interaction potential 
is replaced by the two-body T-matrix can be found in 
24, 31, 56,|6l[7i[7l|73. This replacement naturally 
involves a renormalization of higher order terms. 

Rather than rederive these results, we simply intro- 
duce the T-matrix into the theory from the start, which 
is why the contact potential appears in Eq. ^ for the full 
quantum field operator. This does mean, however, that 
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we have implicitly included at the outset various physical 
effects (multiple two-body scattering) which must also 
appear in the many-body treatment. To avoid double 
counting these effects we must identify where they ap- 
pear in higher order calculations and explicitly subtract 
them in the appropriate approximation. Since the lead- 
ing order contribution from interactions is the nonlinear 
part of the GPE, this is the first term requiring renor- 
malization. In the GGPE, the interaction strength Uq in 
this term must therefore be replaced by the expression 

mill 

Uo^Uo + A[/o, (47) 

where AUq is the second-order approximation to the 
interaction strength as calculated from the Lippmann- 
Schwinger equation 

/d^k 1 

We note that this correction is appropriate even for 
gases in a trapping potential despite the fact that it in- 
volves an integral in momentum space. This is because 
Uo corresponds to the homogeneous T-matrix, so for a 
trapped gas two-body collisions have been implicitly in- 
cluded via a local density approximation and must there- 
fore be subtracted in the same manner. The integral in 
Eq. (|48() is itself UV divergent, but this cancels with the 
divergence in the anomalous average so that the combi- 
nation is finite [see Eq. (|^ ]. 

We therefore have a correction term in the GGPE 
which is iVoA[/o|*i'2(r, i)P<l'2(i", i)- This correction can 
be grouped with the diagonal part of the anomalous av- 
erage, leading to a finite, renormalized version of this 
quantity defined by 

m«(r, t) = m(r, t) + No^^^r, t). (49) 

The renormalized, generalized GPE is therefore given by 
Eq. ifTUIl with m(r,t) replaced by m^(r, i). 

We note that this procedure only affects the diagonal 
part of the anomalous average. The off-diagonal terms 
TO(r, r', t) with r 7^ r' do not require any UV renormaliza- 
tion and are unaffected. In particular, the contribution 
to /(r, t) in Eq. from ?7i(r, r', t) does not contain any 
UV renormalization. 



IV. THEORY OF LINEAR RESPONSE AT ZERO 
AND FINITE TEMPERATURE 

In this section we apply the number-conserving, mean- 
field theory described above to consider the linear re- 
sponse of a Bose-condensed system at zero and finite 
temperature. We closely follow the method of Giorgini 
m, although the calculation differs in detail because it 
contains effects due to the orthogonal projector Q and 



the term /(r, t) in Eq. Hll() which are absent in the more 
usual symmetry-breaking formalism. We have found nu- 
merically that these terms can have a significant effect 
for the parameters of the JILA experiment 0, |^ |6| . We 
also explicitly include the effect of the external perturba- 
tion on the non-condensate. This is generally neglected in 
analytical calculations but is vital for an understanding 
of the JILA experiment 0j as well as for an accurate 
description of the dipole modes. 

We start in Sec. IIV Al with the simplest theory ob- 
tained by using only the GPE, which is generally suffi- 
cient at zero temperature. In Sec. IIV 51 we then consider 
the corrections which arise from the GGPE. These are 
important at finite temperature but they also lead to 
changes in the zero temperature results. Our results ap- 
ply to the collisionless limit in which the characteristic 
time scale for quasiparticle collisions is long compared 
with the chemical potential. More formally, the theory 
is based on a systematic expansion in a small parameter 
which for a homogeneous system of volume V is given by 

M 

(na3)i/2 = 0), 




where no — No/V and n — N/V are the condensate and 
total densities respectively. The origin of these expres- 
sions is discussed further in Sees . IV Al and IV B1 

We note that the theory is not an expansion to second 
order in the coupling constant Uo as stated in 28] . This 
can be seen from the fact that terms of order Uq, cor- 
responding to collisions within the non-condensate, are 
neglected. Instead, we calculate the leading order cor- 
rections to the Bogoliubov theory in the small parameter 
given above. We show in Sec. IV Bl that non-condensate 
collisions are of higher order in this parameter and can 
therefore be neglected for the purposes of this calculation. 
A kinetic theory which includes all terms of order C/g has 
been given by Walser et al. 15 4l and, in the pcrturbative 
limit, by Rusch and Burnett [31[. 

We note that at finite temperature the small param- 
eter of the theory is proportional to l/y/no and hence 
becomes large close to the BEG phase transition. Set- 
ting this parameter equal to one gives an estimate for the 
boundary of the critical region which is closely related to 
the Ginzburg criterion for the failure of mean-field the- 
ory |66j |. Inside the critical region perturbation theory 
will fail, but for a dilute gas this region is very narrow 
and our approach is valid for a wide temperature range. 
If Eq. H50() is used to estimate the validity of the theory 
for a trapped gas by replacing no with no(r) via a lo- 
cal density approximation, the small parameter becomes 
spatially dependent and diverges at the edge of the con- 
densate. This is primarily a consequence of the crudity of 
the estimate but perhaps indicates that the theory may 
be more successful for modes localized near the centre 
of the condensate than around its edge. Evaluating the 
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parameter using the central condensate density probably 
provides a more useful diagnostic, but the validity of the 
theory in this case is best checked a posteriori simply by 
seeing if the predicted energy shifts and widths are small 
compared to the unperturbed values. 



A. GPE theory 

The excitations of a Bose gas can by studied by consid- 
ering the linear response of the system to an infinitesimal 
time-dependent external perturbation P(r, t) 25, 26, 28'|. 
The leading order theory (Bogoliubov theory) is obtained 
using just the GPE of Eq. 0. We consider the case that 
the external potential V{r,t) can be written as a sum 
of two terms, a static trap potential Vtrap(r) (if any) and 
the perturbation P(r, t), both of which are real functions. 
We therefore write the GPE of Eq. ® as 

i;i^$o(r,t)= [^sp(r)+P(r,i)-AoW] $o(r,t) 

+iVoC/o|$o(r,t)|2$o(r,t), (51) 

where /fsp(r) is given in Eq. H24|) and A^o is independent 
of time. 

In the static case where P{r, t) = 0, the GPE has time- 
independent solutions <l>o(r, t) — ^o{r) given by Eq. (|23|l . 
We will consider the situation where a condensate has 
been formed at low temperature and has settled into a 
solution of this equation. We note that the formalism we 
will develop applies in principle to any such solution and 
not just the ground state. In particular the results can 
be applied to vortex states, and for this reason we have 
been careful not to assume that $o(r) is real. 

When the external perturbation P(r, t) is applied, the 
system responds with a time-dependent oscillation of the 
mean field 

$o(r,t) = $o(r)+<5$o(r,i). (52) 

In general 5$o(r, t) will include a contribution which rep- 
resents a global, time-dependent phase applied to <i>o(r). 
We retain the freedom to choose this phase by including 
a variation in the parameter Ao {t) 

Xo{t) = Xo+6Xo{t). (53) 

Substituting these expressions into Eq. (|51|l and lineariz- 
ing with respect to (5<l'o(r, t), SXo{t) and P{r, t) gives the 
equation of motion for the fluctuation 

ih^^5^o{v,t) = 4°P)(r)5<fo(r,t) + Mf P'(r)^$*(r,t) 



+ [P(r,t)-ao(i)]$o(r), 



(54) 



where the quantities V"^^\v) and Mj^^\v) are defined 
simflarly to Lo(r) and Mq{y) of Eqs. (EHJ and ^F7^ but 
with the orthogonal projector replaced by unity 



4'^''^(r) = i/,p(r) - Ao + 2iVo^7o|$o(r)|^ (55) 
MfP)(r)=7Vof/o$g(r). (56) 



Combining Eq. (|54() with its complex conjugate and 
taking the Fourier transform, we obtain 

-[P{v,u)-8Xo(u;)\Xe,{v) = Q, (57) 
where the condensate fluctuation spinor I?"I>o(r, w) is 



5$o(r,cj) 
<5$5(r,-c.)/ ' 



(58) 



and C\^^^ (r) is defined analogously to Cq (r) 



(GP) 



(r) 



f (GP) 



(r) 



MfP)(r)' 

^(GP), 



(59) 



The Fourier transforms are defined by 6^o{r,t) = 
c?w 5$o(r, w)e~*'^*, and we have used the fact that 
the components of d^o{r,t) and 6^Q{r,t) oscillating at 
the frequency u; are 6^o{r,uj) and 6^Q{r, —uj) respec- 
tively. 

We solve Eq. H57|l by expanding the fluctuation in the 
complete set of states provided by the eigenstates of Co{r) 
given in Eq. 



(60) 



where the sum is over all states (i.e. with both pos- 
itive and negative norms, zero and non-zero energy). 
The form of V^o (r, (^) and the relation between posi- 
tive and negative energy spinors of Eq. H31|l implies that 
&_j(cj) = b*{-uj) [in the time-domain 6_i(t) = b*{t)]. 
The expansion coefficients {bi^Lu)} can be projected out 
of this expression using the orthonormality relation of 
Eq. H30() . Specifically, we have 



drA;t(r)a3p$o(r,w). 



(61) 



It is convenient to separate the contribution from the 
two zero-energy condensate modes, and write the expan- 
sion in terms of the norm and phase states defined in 
Eq. H35() . Hence we write 

I?$o(r, w) = bN„iuj)XN„{r) + ibe„iuj)Xeoir) 

+ ^ 6,(^)A',(r), (62) 

where the sum is now over all non-condensate modes and 
bNo{Lo) and bgg(u}) are respectively the Fourier transforms 
of the real and imaginary parts of b+o (t) (the coefficient of 
the positive norm, condensate mode in the time-domain). 
The first term in this expansion corresponds to changes in 
the norm of the condensate while the second corresponds 
to changes in phase. Conservation of the norm of the 
condensate wavefunction gives 67Vo(<^) = and we will 
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choose the parameter d\o{uj) so that there is no phase 
mode, i.e. bgg{uj) = (see below). 

The action of £Q^^''(r) on the Xi{r) spinors is easily 
found using the fact that these are eigenstates of Lq (r) . 
For 7^ we have 



(63) 



where the coefficients Ci are the time-independent ver- 
sions of Eq. gill with $2(r,i) ^ $o(r) 111 



c. = ydr7Vof/o|$o(r)P[$S(i-)M^(r) +$o(r)«^(r)]. (64) 
For the phase and norm modes we have 



,(GP) 



{r)Xo,iv)^0, 



(65) 



4'^P)(r)A'Ar„(r) = 2iVoC/o|*o(r)pA'e„(r). (66) 

Equation H64|l shows that the a coefficients are only non- 
zero for modes with the same symmetry as the conden- 
sate. In particular, they vanish in the homogeneous limit 
and are only non-zero for gases in anisotropic, harmonic 
traps for modes with positive z-parity and zero axial an- 
gular momentum (m = 0). If Ci = 0, Eq. shows 
that the corresponding spinor Xi{r) is an eigenstate of 

c'^^\r) as well as of £o(r). 

Substituting the expansion of Eq. I|62|) into Eq. I|57l) 
and projecting out the coefficients using Eqs. H61|) and 
gives 



bp{uj) 



(67) 



ihwheaiuj) = Pm{oj) - 5\q{lo) + ^ bi{ijj)ci, (68) 

where the excitation matrix elements P^oioj) and Ppo(^) 
are defined by 

Poo(c^)- JdvP{r,Lu)\<i>o\^ (69) 
Ppoiuj) = jdvP{v,uj) [u;^o + v;^l\ , (70) 

(Poo('^) is just a special case of Ppo('^) corresponding to 
p = ±0). Eq. (|68|l shows that the choice 



5\q{lj) = Poo(^^) + X! 



(71) 



gives bgg{uj) = and so removes global phase fluctua- 
tions of the condensate. The fact that the coefficients 
Ci are only non-zero for modes with the same symmetry 
as the condensate, means that a limited subset of modes 
contribute to (5Ao(w). In the homogeneous limit all these 
coefficients are zero and therefore the issue of condensate 
phase fluctuations does not arise if Poo('^) = 0- 



The coefficients {bi{uj)} describe the contribution to 
the condensate density response from the various quasi- 
particle modes. The condensate density fluctuations 
6nc{r,t) — 6{No\^o{r,t)\'^) can be written in terms of 
these quantities using 

5n,{r, Lu)^NoY, U^) [^^oi^Ki^) + M^Mr)] ■ (72) 

i=i±0 

This expression provides the link between theory and ex- 
periment because experiments generally measure conden- 
sate density fluctuations, whereas the theory gives ex- 
pressions for the coefficients bi{uj). 

The coefficient bp{Lo) in Eq. H67|) diverges if the con- 
densate is driven at the resonance frequency of the cor- 
responding excitation, lo — tp/h. This is, of course, be- 
cause there is no damping in the theory at this level of 
approximation. Finite quantities can be obtained using 
the standard prescription of inserting a small imaginary 
part in the denominator via uj ^ w -\- ij. This procedure 
is essential for numerical calculations at the next level of 
approximation, and can be justified from the finite exper- 
imental observation time, as we discuss in Appendix O 



B. Generalized GPE theory 

For a given experimental configuration (type of atom, 
trapping potential etc) the analysis given above has only 
one parameter, namely the condensate population Nq. 
It is therefore effectively a zero-temperature treatment, 
making no explicit mention of the non-condensate. In 
this section, we extend the analysis to finite tempera- 
ture using the GGPE to describe the dynamic coupling 
between the condensate and non-condensate. The previ- 
ous results acquire an implicit temperature dependence 
because we use basis quasiparticle wavefunctions Xi{r) 
and energies corresponding to the relevant equilibrium 
condensate population, A^o — No{T). However, in the 
following discussion, we do not write the temperature 
dependence explicitly and simply denote the equilibrium 
condensate population by Nq. As discussed at the start 
of Sec. IIVI the theory developed in this section contains 
the leading order corrections to the Bogoliubov theory in 
the small parameter of Eq. H50() . 

The starting point for the higher order calculation is 
the GGPE of Eq. ^ with Hspir,t) = i?sp(r) + F(r,t) 
and UV renormalization as in Eq. 149|) . In the static case, 
the time-independent solutions <I'2(r,t) ~ *I'2(r) satisfy 



^sp(r) - A2 + (A^o + AAo)[/o|$2(r)pJ $2(r) (73) 
+2[/on(r)$2(r) + C/om«(r)$;(r) - /(r) = 0, 



where A2 is a constant, again roughly equal to the chemi- 
cal potential in an improved approximation. n(r), m^(r) 
and /(r) are calculated from Eqs. and us- 

ing static quasiparticle wavefunctions and energies ob- 
tained by solving Eq. ^ with A^q = No{T). 
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As before, application of the external perturbation 
P{r,t) leads to all quantities developing a small time- 
dependent oscillation around their static values. We 
therefore write 

$2(r,<) = $2(r) +(5$2(r,i), 
A2(t) = A2 + (5A2(i), 

N„{t)^No + SNo{t), \ (74) 

n{r, r', t) — h{r, r') + 5n{r, r', t), 
rn^(r,r',<) = m^(r,r') + (5m^(r, r', i), 
/(r,<) = /(r) + 5/(r,t). 

We neglect here any dynamic change in the variance of 
the condensate population AiVg, although changes in the 
mean are included via (5iVo(t). Substituting these expres- 
sions into Eq. H1U|) and linearizing with respect to all 
small quantities, leads to an equation of motion for the 
condensate fluctuation (5$2(r, i)- As before we take the 
Fourier transform of this and combine it with its complex 
conjugate in the two-component form I?$2(r,ti') defined 
in the same manner as Eq. H58I) . 

Linearizing with respect to the quantities in Eq. H74|) is 
exact in the limit P(r, t) — > 0. However, we now perform 
a second linearization which is governed by the small pa- 
rameter of Eq. H50() 



Here ALj^^"* (r) and AMj^^"* (r) are defined by 



*2(r) = $o(r) + A$2o(r), 
A2 = Ao + AA20. 



(75) 
(76) 



A$2o(r) and AA20 are static quantities which arise from 
the difference in the condensate shape and energy be- 
tween the generalized and ordinary GPE descriptions. In 
all terms which were not present in the GPE theory of the 
previous section, we make the leading order replacement 
^2(r) — > <l'o(r) and A2 — > Aq, while in terms which were 
present in the earlier treatment we keep A$2o(r) and 
AA20 to linear order. We note that because these correc- 
tions are independent of P(r, t) we must retain quantities 
such as A$2o(r)(5<I'2(r, t). These are not quadratic in the 
perturbation, but instead are linear in the perturbation 
and also linear in the small parameter of the theory. 

Carrying out this second linearization, leads to the fol- 
lowing equation for the condensate fluctuations 

[hLu - 4''^^(r)] V^2{r,Lo) - [Pir,Lo) - SXl,{Lo)]Xg„{r) 
( A4^P)(r) AA4^P)(r)\ ^ , , 



Sf{r,u;) 



+ SNo{Lu)Uo\Mr)\^Xe,ir). 



AL 



^^^\r) = 2Uon{v) - AA20 + 2AiVoC/o|$o(r)P (78) 
+ 2iVoC/o[*o(r)A$;o(r) + $S(r)A$2o(r)] , 



+27Vot/o$o(r)A$2o(r), 



(79) 



and 5\2{lo) and 5X\{u}) are respectively the Fourier trans- 
forms of the real and imaginary parts of 5\2{t). In fact, 
the imaginary part of 6\2{t) makes no contribution to 
the evolution of a mode with finite energy because it ap- 
pears multiplied by the condensate norm mode which is 
orthogonal to all such states. 

We solve Eq. (|77|) as before by expanding the fluctua- 
tion in the static quasiparticle basis of Eq. (|28|l 

I?<I>2(r,w) = bNa{ijj)XN„{Y) + ihe„{uj)Xe„{Y) 

+ h{uo)X,{Y). (80) 

Unlike the previous expansion of Eq. (|62|) . however, the 
coefficient 6jVo(w) is not zero, even though the norm of 
$2(r,i) is conserved by the GGPE. The reason is that 
the condensate wavefunction appearing in the condensate 
modes of the expansion is $0(1") and not $2(r)- In fact, 
conservation of the norm of $2(1", i) gives 



^JVo (w) 



2 ^ 



bi{uo) j dr (A$;oM, + A$2oWi 



(81) 



(77) 



In this equation we have used the fact that bNg^u) is zero 
at leading order and also that the relative phase of $2 (r) 
and $0 (r) can be chosen so that the phase mode gives no 
contribution. As in the GPE calculation, we can always 
choose (5A2(w) so that beg{Lu) — 0. At leading order this 
requires choosing 6X2(^0) = ^Ao(w) as given by Eq. lf7T)l 
and this expression is sufficient for this calculation. We 
will assume from here on that this choice for SX2{uj) has 
been made. 

Substituting Eq. (|Sn|l into Eq. (|77|l and projecting out 
the coefficient of mode 'p' ^ ±0 gives 

Xpihu; - ep)bp{cj) - Ppoi^) ^ AP^p{u;) + AP^^^\cj) 
-|-Aii;f)(c.) + A4^)(c.), (82) 

where the various quantities on the right-hand-side are 
defined and discussed below. These new terms corre- 
spond to changes in the excitation matrix element and 
energy of the quasiparticle mode 'p' arising from the 
presence of the non-condensate. The changes are of two 
types, static and dynamic, denoted by the superscripts 
(S) and (D) respectively, and correspond to the different 

roles of the thermal cloud. The static term AEp^\uj) 
comes from interactions between a condensate fluctu- 
ation and the equilibrium non-condensate mean-fields 
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n(r), m^(r) and /(r). AP^g^a;) describes the effect of 
changes in the static condensate shape [<i^o(r) *i'2(r)] 
on the excitation matrix element Ppo{uj). Dynamic terms 
describe the effect of the couphng between condensate 
and non-condensate fluctuations. 

Non-condensate dynamics can occur via two mecha- 
nisms; either the non-condensate is driven directly by 
the perturbation or it is driven indirectly via the conden- 
sate. The second possibility is the one usually considered 
in analytical discussions of condensate excitations. This 
mechanism gives rise to the dynamic term IS.Ep^'^ [lj) de- 
scribing the generation of non-condensate fluctuations by 
the condensate and their subsequent back action which 
damps and shifts the condensate excitations. The in- 
clusion of this contribution is required to cancel infra- 
red divergences in the static terms and to obtain a gap- 
less excitation spectrum consistent with the Hugenholtz- 
Pines theorem pj] |. This is demonstrated explicitly in 



|2l|2l|2i| and is discussed further in Sec. lvn 

However, the non-condensate can also be excited di- 
rectly by the external perturbation. The resulting fluc- 
tuations can then couple to the condensate and act, in ef- 
fect, as an additional external perturbation. This process 
therefore changes the excitation matrix element PpQ{uj) 

and is described by the term APp^'' (uj). Direct excitation 
of the non-condensate is generally neglected in analytical 
discussions but turns out to be of crucial importance in 
explaining the results from the 1997 JILA experiment, 
as demonstrated recently 0, Q • It is also important to 
include this effect for consistency with the dipole mode, 
see Sec. and 

The static and dynamic shifts can be written as an 
expansion over the various excited quasiparticle modes 
as 



(83) 



Ai?f)(c.)= ^ AEl^HcuMco), (84) 

( S) 

where the quantities AEpq are independent of frequency 

while the quantities AEp^\uj) have a non-trivial fre- 
quency dependence. Equation H82|) can therefore be writ- 
ten as a matrix eigenvalue equation for non-condensate 
modes {p, q ^ ±0) 



^ Hpq{uj)bg{uj) :^Vpo{uj), 
q=i±0 



(85a) 



Hpgiuj) =xp{f^- ep) 5pg - AElf - AE^J^\lu), (85b) 
Ppo (^) = Ppo {uj) + AP^I^ {u) + AP^^^ {u) , (85c) 

where 5pq is the Kronecker delta. 

This equation can be solved by inversion of the matrix 
H(a;). We simplify the calculation considerably, however, 
by assuming that the perturbation P(r, t) has a symme- 
try and frequency width such that only one quasiparticle 



mode is significantly excited. This is appropriate for ex- 
periments which measure the excitation energies of the 
system as a function of temperature. In addition, if we 
are only interested in frequencies near a particular un- 
perturbed resonance to « ep/?i, we can use a rotating- 
wave-approximation and neglect the contribution from 
the mode with energy —ep. In this case, we can as- 
sume that only a single coefficient bp{uj) is non-zero and 
Eq. (|85() can be solved straightforwardly. Introducing the 
experimental resolution via w — *■ w -I- «7 as discussed in 
Appendix |21 we find that the Fourier transform of the 
response coefficient can be written as 



bp{ijj) = XpPpa{oj)np{uj + 17), 



(86) 



where the generalized response function TZp{uj) is defined 
by 



TZp{uj) = 



AP^',\u;) + AP^^\co) 



pO 



Ppo(t^) 



Gpiiu), (87) 



and the resolvent Gp{u!) can be written in terms of a 
frequency-dependent self-energy as 

1 



hu! — ep — 2^p[UJ) 
i:p{u;)=Xp\AElP +AEl^H^) 



(89) 



C/p(ijj) describes the effect of the static and dynamic 
coupling between the condensate and non-condensate in 
the limit that the external perturbation only excites the 
condensate. The additional factor in TZp{uj) accounts for 
the change in the excitation matrix element when non- 
condensate fluctuations are excited directly by the exter- 
nal perturbation. The self-energy 'Ep{uj + ij) is a complex 
function, the real part describing a shift in the energy of 
the excitation and the imaginary part giving its damping. 

If Sp (lu + ij) and AP^^^ {uj + ij) are roughly inde- 
pendent of frequency, the energy shift can be calculated 
straightforwardly from the poles of Gpitu + 17) by finding 
the solutions to 



Ep = fkjp = Real[ep + T,p{ujp + 17)] . 
The corresponding decay rate is then given by 



Tp = — Imaj 



[T,p{ujp + ij)]/fi. 



(90) 



(91) 



This situation arises when an excitation couples to a con- 
tinuum of decay channels, as in the homogeneous limit, 
and the resolvents Gp and TZp are Lorentzians. For a finite 
system, however, I]p(w -I- 17) depends on frequency, and 
neither Gp nor TZp are perfect Lorentzians. In this case 
the line shape depends on the details of the system, in 
particular the availability of decay channels, the coupling 
strength to each of them, and the broadening of the levels 
involved. We can still extract energies and decay rates, 
however, by fitting hp{uj) to a complex Lorentzian (the 
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experimental resolution 7 should be subtracted from the 
resulting decay rate). An alternative is to mimic the ex- 
perimental procedure exactly and Fourier transform ap- 
propriate moments of the condensate density response 
(Eq. below) to the time-domain and fit a decaying 
sinusoid. 

As before, the coefficients {bp{uj)} can be used to con- 
struct the condensate density fluctuations, which are the 
experimentally relevant quantity. At this order, how- 
ever, we must include both the dynamic fluctuations in 
the condensate number SNo(lu) and also the non-zero co- 
efficient bNg{LLj). Equation H72|l is therefore replaced by 

Sn,{r,u) = [SNoiu) + 2NobNo{L^)] |*o(r)P (92) 
+ NoY^ h{uj) [%{t)u,{v) + $o(r)w,(r)] , 

where 5Nq{lo) is defined by 5Nq{uj) = — Jdr(5n(r, w) and 
is given explicitly in Eq. I|A38|I . We can also calculate the 
spatial dependence of the non-condensate density fluc- 
tuations 5n{r,uj) using the explicit expression given in 
Eq. ljA26|) . This can then be used to investigate the issue 
of the relative phase of condensate and non-condensate 
oscillations [Ulillil. 

We now define the various quantities appearing in 
Eqs. (|82() - (|84|l and discuss their physical interpretation. 
Details concerning the derivation of the expressions are 
given in Appendix IXI 



1. Definition of static shifts 

We define static shifts as those terms which do not 
involve any dynamic non-condensate fluctuations. The 
only frequency dependence in these terms therefore arises 
from the expansion coefficients {bi{Lo)}. 

The energy shift AEpq'' of Eq. H83|) can be written as 



= AE^ip, q) + AExip, q) + A£;shapo(p, q) 

+ AEANoiP,q) + ^Ef\p,q), (93) 

where the various contributions are defined by 

AEiip, q) = jdr ^2Uoh + v^Vq] 

+Uom^ulvq + Uom^*vluqY (94) 



AEx{p, q) = -AA20 Jdr [u*pUq + v*Vq] , 



A-Eshapo(p, 9) = 

2N0U0 jdv[ [$oA$;o + *o A3'2o] [u;uq + v;vq] 

- jdr[cl[A^l^Uq + A^2^Vq] 

+ c,[A$2o^/; + A$>;]}, (96) 

AEANoiP,q) ^ ^NoUoJdr ^^2\<i>of + 

*0 VpUqj, 



'^lu*pVq 



(97) 



AEf\p,q) 



^ (A,.o +^2,o.) ^ ^ ^ 



j>0 



E 

i>0 



2VNo 

{Blqio + B2qio) 



[N^Up, + (N, + l)Vp, 



(98) 



( 91 

The various quantities in AEj {p, q) are defined by 

U,j= [dru*{r)u,ir), (99) 

(100) 



Vij = drv*{r)vj{r), 



Wij = I drui{r)vj{r), 



(101) 



(95) 



while Aqio, Blqio, B2qio and B2qOi are special cases of 
the coefficients Aqij, Blqij and B2qij defined later in 
Eqs. (|1L)7|I - I|109|) [the index '0' refers to the condensate 
mode with positive norm of Eq. (|32|) ]. 

The quantities AA20 and A$2o(r) needed to evalu- 
ate AE\{p,q) and AEshapc{p,q) are calculated in Ap- 
pendix ^3 Briefly, the static solution to the GGPE of 
Eq. H73|l is obtained by linearizing around the static so- 
lution to the ordinary GPE using Eqs. I|75() and l|76|) . 
The resulting equation is solved using an expansion in 
the quasiparticle basis of Eq. H28|) . This linearization is 
necessary for consistency with the linearized treatment 
of condensate shape and energy effects in Eq. lf77|) and 
also for a systematic treatment to a given order in the 
small parameter of the theory. 

The static change in the excitation matrix element 
AP^^\uj) is defined by 

AP^piL,) ^ Jdr [P(r,c.)-Poo(^)] KA$2o + A$y . 

(102) 

The contribution from Poq{uj) arises from the choice of 
6X2(10) = SXq{ll!) given in Eq. H71|) which removes con- 
densate phase fluctuations to leading order. The rest of 
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the expression is a straightforward modification of Ppo (uj) 
as given in Eq. H7U|) accounting for the change in shape 
of the condensate. 



2. Interpretation of static shifts 

We now discuss the physical interpretation of the vari- 
ous terms appearing in the static shifts. Ai?4(p, q) is the 
energy shift arising from the interaction between the con- 
densate fluctuation and the static mean-fields of the non- 
condensate n(r) and TO^(r). In an alternative approach 
to the theory of the Bose gas, based on perturbation 
theory with the many-body Hamiltonian, this contribu- 
tion arises from the use of first order perturbation theory 
on products of four non-condensate operators (hence the 
choice of notation) |2Jj . Aijj-'^-' (p, q) arises from the defi- 
nition of the non-condensate as orthogonal to the conden- 
sate while AEano {p^ q) expresses the effect of the vari- 
ance in the condensate population on the nonlinear term 
in the GPE. 

AEx{p,q) and A£'shapc(p, ?) both arise from changes 
in the properties of the equilibrium condensate wavefunc- 
tion due to interactions with the static non-condensate 
mean-fields. The condensate eigenvalue A is the zero of 
energy for excitations so a change in this gives rise to 
AE\{p,q). The first two lines of AEshapdp, q) have a 
straightforward interpretation, describing the change in 
the condensate mean-field potential when the condensate 
shape changes. The final two lines come from a more 
subtle effect. In the current formulation these arise from 
the non-zero condensate norm coefficient b^^luj) given in 
Eq. (|81|l and also from the choice of 5X2(1^) which removes 
pha se fluctuations (see Appendix lA 21 for the details). In 
|24| , the same contribution is obtained from the fact that 
quasiparticle excitations are defined to be orthogonal to 
the condensate. A change in the shape of the conden- 
sate therefore leads to a change in the definition of the 
orthogonal subspace and this in turn affects the energy 
of the excitations. 



3. Definition of dynamic shifts 

We define dynamic terms as those which involve fluc- 
tuations in the non-condensate. These terms have a non- 
trivial frequency dependence over and above that con- 
tained in the expansion coefficients {6^(0;)}, as indicated 
in Eq. (jH3. 

Expressions for the energy shift AEp^\oj) and the 
change in the excitation matrix element AP^^"' (lu) are 
derived in Appendix lA .31 Briefly, the derivation con- 
sists of writing h{r,r',t) and m(r, r',t) in terms of the 
time-dependent Bogoliubov functions Xi (r, t) and evolv- 
ing these according to Eq. (|40|l . Linearizing the changes 
in these quantities around their static values and solving 
the resulting equations in terms of the static Bogoliubov 



basis, leads eventually to expressions for Sh{r,r',Lu) and 
Srhir^ r', uj) from ^ 
straightforwardly. 



Srhir.r' ,uj) from which AEp^\uj) and APp^\uj) follow 



The resuh for AE^pg\uj) is 



y{A)*y{A) 

[1 + N, + N,] 



2{fuj + e^ + ej) 



ij>0 

E 

ij>0 



y(Bl)*y(Bl) 

pij qij 
2 {huj — ei — tj) 



[1 + N, + Nj 



yiB2)*y(B2) 

pij qij 



+ AEo{p,q) + AE''ip,q), 
where AEq (p, q) is defined by 
^Eo{p,q) = 

--^ E [Y^of '*W.qN. + Yi^^*W.,{N. + 1) 



(103) 



i>0 



'=Y.[Y:;T*u.,N.+Yir*v.,iN.+i) 

° i>0 



(104) 

and AE^{p, q) is the UV renormalization 
AE^ip, q) = 2N0AU0 Jdr |$o|' + , (105) 

with AUq as in Eq. (jlSll . 

The coefficients Yp^\ ^'^'^ ^if^^ ^'-^ defined by 



^p J 



/iVo 
c* 



J 7 - 



Jr) 



' Jp^, (106a) 



/No 

^j J* 



(106b) 



(106c) 



where the coefficients Ci are given in Eg. I[64ll and Apij, 
Blpij, B2pij, lij and Jy are defined by |79| 

Apij = 2y^UQ Jdrup[^l {uiVj + ViUj) + ^oViVj] 



+Vp[^0 {UiVj + ViUj) + ^gUiUj] 



(107) 



Bl 



pv 



2^/n^Uo fdvup[^l {u,Vj + ViUj)* + $0^ 



+vp [$o {uiVj + v^u^)* + <PoV*v*] , 

(108) 
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+Vp [$o «Uj + V*Vj) + %V*Uj] , 

(109) 



hj = I dr [u*Uj + v*Vj 



Jij = Jdr [u.iVj + v^Uj] = 2Wt. 



(110) 
(111) 



The term AEoip,q) in Eqs. ifTH!^ and Cnijl, is the 
exphcit contribution to AEpj^\uj) from the condensate 
modes (note that the summations in Eq. (|103() exclude 
the state +0). These modes appear in the description of 
6n{r,r',uj) and 6m{r,r',uj) because the time-dependent 
Bogoliubov functions Xi{r, t) are orthogonal to the time- 
dependent condensate wavefunction (see Appendix I A 3p . 

The dynamic change in the excitation matrix element 
APp^' (w) has a form similar to AEp^^ {lj) 



ij> 



2y/No {hiu + + e^) 



[1 + N, + N^] 



y(BI)* p(B) 



ij>0 



2y/No (hio - - ej) 



y(B2)*p(L), . 

E [iv, - AT,] , (112) 



y>0 



'/Vo(fiw - ei + ej) 



where we have defined Beliaev and Landau perturbation 
matrix elements by 

P^piu;) ^ Jdr [P(r, uj) - Poo(c^)] [u*v* + v^u*] , 

(113) 

P,f \co) = [dv [P{r,u;) - Poo(^)] [u*u,+v*v,], (114) 



with Poo{uj) as in Eq. (|69|l . To the best of our knowl- 
edge Eq. (|112|l has not been obtained previously in the 
literature. 

To calculate the condensate density fluctuations in 
Eq. H92|l we also need an expression for SNo{uj). This 

has the same structure as AEp^\oj) and APpl^\uj) and 



is written explicitly in Eq. HA38I 



4- Interpretation of dynamic terms 

The self-energy AEp^\uj) describes the effect of the 
dynamic coupling between the condensate and non- 
condensate in the coUisionless limit. In this approx- 
imation, the non-condensate is modelled as a weakly- 
interacting gas of quasiparticles and its dynamics there- 
fore corresponds to rearrangements of existing quasipar- 
ticles among the available levels (Landau processes) or 



the creation or annihilation of quasiparticles (Beliaev 
processes). More specifically, the first line of Eq. H103|) 
corresponds to the collision and annihilation of three 
quasiparticles, while the second line describes the spon- 
taneous decay of a single quasiparticle into two others. 
The third line corresponds to Landau processes in which 
two quasiparticles collide and coalesce to form a third. 

The matrix elements Y^f^ ^pij^^ ^^'^ ^pij^^ describe the 
coupling strengths for these processes and contain the rel- 
evant selection rules (such as angular momentum conser- 
vation and parity) . If the quasiparticle modes are ther- 
mally populated these processes are enhanced by Bose 
stimulation of the various scattering events involved. The 
dynamic shift thus corresponds physically to the intro- 
duction of many-body T-matrix effects into condensate- 
non-condensate collisions [the many-body T-matrix is 
introduced into condensate-condensate collisions by the 
static anomalous average rh^{r)] [M HI IM El IzS Izl ■ 
Landau and Beliaev processes can occur in a virtual 
sense or in a real sense when one of the energy denomina- 
tors vanishes. To obtain finite quantities in calculations, 
it is therefore essential to include a small imaginary part 
in the frequency via uj ~* uj + ij, as in Eq. (|86l) . This is 
justified by the finite experimental resolution as shown in 
Appendix[ni AEpj^\uj + 17) is then a complex quantity, 
the real part describing the shift in the energy of a con- 
densate fluctuation and the imaginary part describing its 
damping. 

We now consider the structure of AeI^\ll>) in more 
detail. The numerators in Eq. H103|) contain a product of 
two of the Y matrix elements. The various contributions 
to these terms, given in Eq. H106|) . have different physical 
origins depending on whether they appear with the label 
'pij' or 'gij'. For example, the term involving Cp in the 
Ypij factors comes from condensate number fluctuations, 
i.e. from 5Nq{u}) in Eq. |(77J). The corresponding term 
involving Cq in the Yqij^s instead comes from condensate 
phase fluctuations as shown in Appendix lA 31 There is 
thus an elegant symmetry in the way the theory deals 
with condensate phase and number fluctuations. 

Similarly, the contribution from the terms involving Ci 
and Cj in the Ypij factors comes from the dynamic part 
of Sf{r,uj) in Eq. ||77J), while the corresponding terms 
in the Yqij factors come from changes in the orthogonal 
projectors Q in the equation of motion for the excited 
states. This demonstrates the connection between /(r, t) 
and the issue of orthogonality between the condensate 
and non-condensate. The matrix elements Aqij etc in 
the Yqij^s describe the driving of the non-condensate by a 
condensate fluctuation in mode 'q', while the correspond- 
ing terms in the Ypij factors describe the back action of 
the non-condensate on the evolution of mode 'p'. 

Although it is physically obvious that condensate num- 
ber fluctuations can change the energy of an excitation 
it is perhaps less clear that phase fluctuations can have 
a similar effect. In this regard we note that a global 
change in the phase of both the condensate and all the 
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excited states has of course no physical consequences. 
If we put X2{t) — > A2(i) + dX where dX is a constant, 
then it is easy to see from Eqs. (tTHIl . Ip!?^ and lB?Hl 
that <i>2(r, i), Ui{r,t) and Ui(r, i) transform according to 
$2(r,i) ^ $2(r,i)e*''^*/'', u^{r,t) u,{r,t)e"^^^^^ and 
Wj(r,t) Vi{v,t)e^'^'^^*-/^\ Comparison with Eqs. lfTn|l . 
(|42|l and H43|) shows that the additional phase factor 
cancels in the GGPE. However, the perturbation does 
not necessarily cause the same phase fluctuations in the 
condensate and the excited states so there may be rela- 
tive phase fluctuations which are physically meaningful 
and which can affect the energy. For example, the term 
Poo(w) in the Beliaev and Landau perturbation matrix 
elements of Eqs. and H114(l comes from conden- 

sate phase fluctuations [c.f. Eq. H71|l ]. This term is sub- 
tracted from the perturbation driving the excited states 
so that only the relative phase of the condensate and 
non-condensate has any physical effect. 

As shown by Eq. and noted below Eq. (I7TJ, the 
Ci coefhcients are only non-zero for modes with the same 
symmetry as the condensate. This means that conden- 
sate phase and number fluctuations will only have an 
effect on the response if one or both of the modes 'p' and 
'q' has the appropriate symmetry. Similarly, the effects of 
5f{r,uj) and the orthogonal projector Q in the Ypij/Yqij 
coefficients are only relevant for subspaces where either i 
or j has the symmetries of the condensate. These spaces 
form a small subset of the whole, which means that these 
contributions describe finite size effects and should be 
negligible in the thermodynamic limit. Nonetheless, they 
can be significant for the finite systems of recent experi- 
ments. This is discussed further in [l5Ll6fi|| where we con- 
sider the numerical implementation of the current theory. 
All these effects vanish in the homogeneous limit where 
the Ci coefficients are identically zero (see Sec. IV A|) . As 
far as we know, the contribution to the dynamic shift 
from the Cj's has not been obtained previously. 

The quantity APpq\uj) given in Eq. H112(l describes 
the excitation of the non-condensate by the perturba- 
tion and the subsequent coupling to the condensate. The 
structure of AP^q^ (w) is similar to that of AEp^^ (w) and 
the physical processes involved are the same, the only 
difference being how the non-condensate fluctuations are 
created. In this case, the matrix element P^^j (lj) de- 
scribes the creation of pairs of quasiparticles by the per- 
turbation while P^^^ (w) describes the rearrangement of 
existing quasiparticles. The Ypij factors describe the sub- 
sequent coupling to the condensate. However, APp^\uj) 
has no contribution from the condensate modes and has 
no UV renormalization. 

We note that the factor of 1/^/lVo in the definition of 
APpQ^^ {to) in fact cancels with a factor of ^/Nq contained 
in the definition of the Ypij coefficients [see Eqs. (|lU6f) - 

((TTT|l and Eq. lIHlll]. This means that AP^^\uj) has no 
explicit dependence on Nq, although of course it still has 
an implicit dependence on this quantity via the quasi- 
particle energies and wavefunctions from which it is con- 



structed. It is therefore expected to be only weakly- 
dependent on A^o for a given temperature. 



V. DISCUSSION 

A. Homogeneous limit 

The theory presented above simplifies considerably in 
the homogeneous, thermodynamic limit where a num- 
ber of its properties can be studied analytically. In this 
limit, the static solutions to the ordinary and general- 
ized GPE are simply <i>o(r) = <I'2(r) = 1/VV where 
V is the volume of the system, so A<I>2o(r) = and 
there are no condensate shape effects. The static quasi- 
particle solutions of Eq. I|28|) are simply plane waves 
<%k(r) oc e^^ '' with the famous Bogoliubov dispersion re- 
lation e/c = [(cfifc)^ + (/i^fc^/2TO)^]^/^, where the speed 
of sound is = noUo/m. The parameters Ci defined 
in Eq. (|64|l are zero for i ^ ±0 and therefore the Y^'^^ ^ 

Y^f^^ and Y^^^^ coefficients defined in Eq. (|106|l sim- 
plify to the integrals Apij, Blpij and B2pij of Eqs. H1Q7|I - 
(|109|l . It also follows that the equilibrium value of /(r) 
is zero. Using the expression derived by Giorgini and 
co-authors for the variance in the condensate popula- 
tion 1^31, it is straightforward to show that terms in- 
volving AA^Oi vanish in the thermodynamic limit. The 
energy shift AEf\p,q) of Eq. ^ and the contribu- 
tion from the condensate modes to the dynamic shift, 
AEo{p, q) of Eq. (|104|l also vanish in this hmit. Thus the 
additional terms introduced by the number-conserving 
approach are finite size effects and are not relevant in 
the homogeneous, thermodynamic limit. The theory 
therefore reduces in this case to the Beliaev-Popov the- 
ory which has been discussed by a number of authors 

The remaining expressions for the static and dynamic 
energy shifts can be calculated explicitly. The excited 
states form a continuum, so the self-energy Sp(w -I- 17) 
is only weakly dependent on frequency over the width 
of a resonance and excitation energies and decay rates 
can be calculated using Eqs. (|90|l and H91|l and taking 
the limit 7^0"^. The Beliaev and Landau decay rates 
calculated in this way generate all the standard results 
in the literature [H IH El S H [H [s?, "eO^ and are 
proportional to the small parameter of Eq. (|50|l . 

When the energy shifts are calculated it is found that 
the static and dynamic terms are both large, each giv- 
ing an infra-red divergent contribution to the energy of 
the form A£;^^^ ^^E^^^ cx 1/fc as k 0. How- 
ever, as shown explicitly in [23. l28l | . the divergence in 
the static part is exactly cancelled by the divergence in 
the dynamic part and the resulting energy shift is finite, 
gapless and proportional to the parameter of Eq. 150|l . 
The validity condition for the theory then comes from 
requiring that the energy shifts and decay rates are small 
compared with the unperturbed Bogoliubov energies. 
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Non-condensate collisions 



The theory developed in Sec. IIV Bl represents a sys- 
tematic treatment of a dilute Bose gas in the collision- 
less limit when the validity conditions of Eq. H5l)|l are 
satisfied. As the density and temperature increase, how- 
ever, we will eventually reach a regime where higher or- 
der terms become important and ultimately the theory 
will fail in the hydrodynamic limit. The most important 
physical processes we have neglected are collisions of two 
non-condensate atoms where each remains in the ther- 
mal cloud. Of course a rigorous treatment of this process 
would also have to include a variety of other effects of the 
same order |54j . In this section we simply estimate the 
importance of such collisions in the homogeneous limit 
and show that they are of higher order than we consider 
in the small parameter of the theory. 

We can estimate the non-condensate collision rate at 
temperatures greater than the chemical potential using 
classical kinetic theory. In this case, the collision rate 

^ nc 

is given by 



nc ~ nav, 



(115) 



where a = Sna'^ is the scattering cross-section and v = 
(8fcBr/m7r)i/2 is the mean speed. The non-condensate 
density can be approximated by the critical density n « 
n„.it = C(3/2)/Ai^B where AdB = {2TTh'^ /mkBTY^^ is the 
thermal dc Brogiie wavelength and C{x) is the Riemann 
zeta function with C(3/2) ~ 2.612. Combining these re- 
sults gives 



C(3/2)8ma2(fcBr)2 



(116) 



The principal assumption of our treatment of non- 
condensate dynamics is that quasiparticles form a good 
basis, i.e. they evolve mainly under the influence of the 
condensate mean-field. This means that the time-scale 
for quasiparticle collisions must be long compared to their 
frequency at the relevant energy scale. As the dominant 
contribution to the self-energy comes from quasiparticles 
near the chemical potential fi — uqUo, the importance of 
non-condensate collisions is described by the dimension- 
less parameter 



UqUo 



32C(3/2) 



UqUq 



Ka'). 



(117) 



Comparison with Eq. (|50|l shows that this is proportional 
to the square of the small parameter of the theory, prov- 
ing that non-condensate collisions are only relevant at 
the order beyond the current calculation. 



C. The dipole mode 

The dipole mode (also known as the Kohn mode) cor- 
responds to a centre-of-mass oscillation of the atomic 



cloud, with the condensate and non-condensate oscillat- 
ing in phase. In a harmonic trap, this means that the 
atoms see a time-dependent potential which is a linear 
function of position across the cloud, corresponding to a 
spatially uniform force. Every atom therefore undergoes 
the same acceleration so that the relative motion of a 
pair of atoms is entirely unaffected. The evolution in the 
centre-of-mass frame is consequently the same as if the 
cloud was stationary, completely independent of temper- 
ature or particle interactions (this is the essential physical 
content of the generalized Kohn theorem) jsj] . For this 
reason the dipole modes are a very useful experimental 
diagnostic for determining the trap frequencies, but they 
are also an important test of the theoretical description of 
non-condensate dynamics. In this section we show that 
our linear response treatment describes the dipole modes 
correctly to within the order of the calculation. For this 
to be the case, it is important to include the effect of 
the external perturbation on the non-condensate, via the 
static and dynamic shifts APpg^(tj) and APpq\uj). We 
focus on oscillations along the z-axis but similar com- 
ments apply to the dipole modes in the x-y plane. 



1. The dipole mode in the GPE 

As is well-known, the dipole mode is obtained exactly 
in the GPE of Eq. ® . It is easy to show that for 
a harmonic potential, the time-dependent GPE has ex- 
act solutions consisting of eigenstates oscillating along a 
principal axis of the trap with no change of shape |82l| . 
A consequence of this is that the BdG equations for the 
excitations, Eq. (|28|) . have solutions with energies corre- 
sponding to the principal trap frequencies. 

We can generate dipole oscillations by varying the po- 
sition of the centre of the harmonic trap. For oscillations 
along the z-axis we define P(r, t) by 



y(r, t) = V{t) + P{r, t) = V{x, y, z - Zp), 



(118) 



where Zp{t), the position of the trap minimum in the 
z-direction, is an arbitrary function of time. For a har- 
monic trap of frequency ujz in the z-direction, the spatial 
dependence of the perturbation is therefore given by 



P(r,t) 



,(t). 



(119) 



The condensate responds to this perturbation by os- 
cillating along the z-axis with no change in shape. We 
can therefore write the time-dependent condensate wave- 
function as 



(120) 



where $o(r) is an eigenstate solution of the time- 
independent GPE while zo(i) and hk{t) describe the 
centre-of-mass displacement and momentum respec- 
tively. We note that in general the condensate posi- 
tion does not coincide with the centre of the trap, i.e. 
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zo{t) Zp{t). The extra energy associated with the os- 
cillation changes the rate of condensate phase rotation 
and can be dealt with by choosing the parameter Ao(t) 
via Ao(t) Ao +5\o{t). For a harmonic potential, these 
expressions provide a solution to the full time-dependent 
GPE provided that z^it), k{t) and 6Xo{t) satisfy 



d'^ZQ 
S\o = 



n, hk dzQ 

^uj^{zo-Zp), — = —, 



m 



dt 



2m 



(121) 



I 2/2 2\ 



The dipole mode is therefore obtained exactly in the GPE 
theory, regardless of the amplitude of the oscillations. 

If we now consider the linear limit we see that (5Ao ~ 
because it depends quadratically on small quantities. 
This is consistent with the general expression for 5\q 
given in Eq. I|71|) . which is only non-zero for modes with 
the same symmetry as the condensate. We then have 



(5$o(r,t) = ~za{t) 



dz 



ik{t)z<^Q{Y) 



(122) 



Using the GPE it is straightforward to obtain equations 
for d^Q{Y)/dz and zi'o(r) and from these to show that 
Eq. has the exact solutions [s^l 



XD±{r) - — 



/9$o(r) 

V 



dz 
d%ir) 



dz 



T z<i>o{r)/z, 
±2*S(r)/4 



2 

osc 



(123) 



where Zqsc = {h/mujz)^^'^ is the oscillator length unit 
in the ^-direction and No is a normalization constant. 
These are the positive and negative energy dipole modes 
with frequencies uid± — ±^2. We note that the orthogo- 
nal projectors play no role in the BdG equations for these 
modes because of parity considerations. 

The condensate fluctuation spinor corresponding to 
Eq. H122|l can therefore be written as 



I?*o(r,w) 



kzi;)XD+{v) 



ikz,. 



XD-{r). 



(124) 



Solving the equations for zq and k in the Fourier domain 
gives 



za{uj) ± ik{uj)zl^^ = 



^UJ^Zp{uj) 



(125) 



A general oscillation of the condensate centre-of-mass 
therefore consists of a linear superposition of the posi- 
tive and negative dipole modes with coefficients which 
are Lorentzians centred on and weighted by the fre- 
quency spectrum of the perturbation. 



e mode in the generalized GPE 



2. The 



The nonlinear term in the GPE plays no role in 
the dipole mode because interactions are unaffected by 



centre-of-mass translations. This is true regardless of the 
particular functional form of the nonlinearity, provided it 
depends only on internal coordinates of the condensate. 
Denoting a general nonlinear term by VF(r, i), all we re- 
quire for the dipole mode is that this function translates 
like the condensate, i.e that it evolves according to 



W{v,t) = W{x,y,z- zo)e 



ikz 



(126) 



where W{y) indicates the static form. The fact that 
the shape of the oscillating condensate satisfies the time- 
independent GGPE in terms of W{v) is then sufficient to 
remove this term from the equation of motion. 

The existence of an exact dipole mode in the GGPE is 
therefore assured if we can write this equation in terms of 
some function VF(r,f) which obeys Eq. (|126|l . We there- 
fore define W^(r, t) by the last three lines of Eq. (|10|l with 
UV rcnormalization, i.e. 

W{Y,t) = (7Vo + AA^o)C/o|$2(r,i)P$2(r,<) (127) 
+2t/on(r, i)$2(r, t) + Uofh^iv, i)$;(r, t) - /(r, t). 

This obeys Eq. itT^ if n(r, i), m^(r, t) and /(r, t) evolve 
according to 



n{v,t) = n{x,y,z- zq), 
m^(r, t) = m^(x, y,z- zaje^''^^ 
f{Y,t)^ f{x,y,z- za)e'''\ 



(128) 



and these in turn are assured if the Bogoliubov spinors 
evolve as 











fcj A',(x,J/,z-zo)e~"'*/^ (129) 



The static functions A'i(r) are defined by Eq. H28|) and 
depend on <f>o(r), the static solution to the ordinary GPE, 
while the time-dependent functions Xi{v,t) evolve ac- 
cording to Eq. (|40|l . It is easy to show from Eq. 14U|I . 
that if the perturbation has the form of Eq. H119|l . the 
condensate wavefunction is taken to be $o(a;, y, z— zo)e*'^^ 
and zo, k and 5\ satisfy Eq. H121|l . then Eq. H129|) is sat- 
isfied. In this case W{y, t) satisfies Eq. (|126|) and we have 
an exact solution to the time-dependent GGPE 



$2(r,t) = <^2{x,y,z- zo)e' 



(130) 



where $2(r) is an eigenstate solution. We note that it is 
important to include the effect of the perturbation P(r, t) 
on the non-condensate for this argument to work cor- 
rectly. 

This is not quite the end of the matter, however, for 
two reasons. First, we solve the static GGPE by lineariz- 
ing the change in shape and energy around ^q{v) and Aq 
as in Eqs. (|7S1) and {TSJ. This means that the above argu- 
ment only applies to order A4>2o (r) . Second, the proof of 
Eq. H129|) requires that Xi (r, t) is driven by an oscillating 
condensate with the same shape as was used to define the 
static values, i.e. $o(r). This requires £2 (r, i) ofEq. H19|) 
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to be written in terms of 4>o(r, i) rather than $2(1", 
Nonetheless, the use of $2(1", i) here is essential for cal- 
culating general energy shifts correctly, as discussed in 
Appendix ^ For centre-of-mass oscillations, however, 
the difference between the two wavefunctions is propor- 
tional to the change in shape of the static condensate. 
This is a small correction of order the small parameter of 
the theory and the resulting error is therefore of higher 
order than we consider. 

The above discussion indicates that the dipole modes 
are not obtained exactly in the current theory, but the 
error will be beyond the order of the calculation. This 
can be confirmed by considering the linearized equations 
we actually solve, i.e. Eq. ((77|)- If Eq. (|129|) is vahd, 
this equation can be simplified considerably using the 
linearized expression for A<i>2o (r) given in Eq. (EU. The 
result is that the error in the dipole mode is of order 
(A$2o)^i which is indeed beyond the order of the cal- 
culation. Interestingly, we also find that if we neglect 
the difference between ^o{r,t) and $2(r,i) in high or- 
der terms, but retain the distinction in the leading order 
terms [i.e. the first line of Eq. (|77|) ] then the dipole modes 
are obtained exactly. 

This analysis applies to the case that we use a full ba- 
sis expansion to evaluate 5$2(r, t), which requires solving 
the matrix problem of Eq. H85() . In practice, the calcu- 
lation is simplified significantly by assuming that only a 
single quasiparticle mode is excited. This corresponds 
to writing I?$2(r,<j-') = bD^{u!)X]j+{r), with XD+{r) as 
in Eq. H123() . In this case we find that the error in the 
dipole mode is now of order A$20 but is localized around 
oj = —Uz and so should be very small in the frequency 
range of interest, namely w ~ w^. Numerical results con- 
firming the accurate calculation of the dipole modes in 
our formalism will be given elsewhere |l5j . 

In summary, a full implementation of the theory de- 
scribed in this paper will obtain the dipole modes cor- 
rectly to the order of the calculation, with errors of order 
(A$2o)^- If we additionally assume that only a single 
quasiparticle mode is needed to describe the oscillations 
then the error is of order A<I>2o but with a small prefactor 
at the frequencies of interest. In both cases, it is impor- 
tant to include the effect of the external perturbation on 
the non-condensate. 



VI. CONCLUSIONS 

In conclusion, we have presented a theory of the linear 
response of Bose-Einstein condensates to external pertur- 
bations at finite temperature. We use a full quasiparticle 
description of the non-condensate and its dynamic cou- 
pling to the condensate. An important physical process, 
usually neglected in analytical calculations, is the direct 
excitation of the non-condensate by the perturbation and 
its subsequent coupling to the condensate. This process 
is explicitly included in our analysis and was shown in 
to be responsible for the anomalous behaviour of the 



m = excitation observed in the 1997 JILA experiment 
It is also necessary to include this process for a cor- 
rect description of the dipole modes. 

The theory is based on the time-dependent, number- 
conserving equations of Castin and Dum -E^ - This for- 
malism describes a number of finite size corrections com- 
pared with the more usual symmetry-breaking approach, 
and these can be significant in calculations for recent ex- 
periments . The dynamic coupling between condensed 
and non-condensed atoms is treated using perturbation 
theory in a quasiparticle basis. Despite statements to the 
contrary in the literature (TsL fl^ , a perturbative calcu- 
lation is appropriate provided that care is taken to ac- 
count for finite size effects and direct excitation of the 
non-condensate. 

The theory developed here is a systematic, gapless ex- 
tension of the Bogoliubov theory which is valid in the 
collisionless limit of well-defined quasiparticles. It is 
consistent with the generalized Kohn theorem for the 
dipole modes and includes the time-dependent normal 
and anomalous averages, Beliaev and Landau processes, 
and all relevant finite size effects. Details on the numer- 
ical implementation of the theory for trapped gases will 
be given elsewhere [T5| . 
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APPENDIX A: ENERGY SHIFTS 

In this appendix we give some details of the derivation 
of the formulae for the static and dynamic energy shifts 
AiJpg-' and Ai5p^''(a;) and the changes in the excitation 
matrix elements AP^p (lu) and APpq \lu) quoted in the 
main text. We also give expressions for the change in en- 
ergy and shape of the equilibrium condensate wavefunc- 
tion. We use the words 'static' and 'dynamic' for energy 
shifts to refer to the role of the non-condensate. In the 
linearization procedure, static shifts therefore arise from 
terms involving the condensate fluctuation (54'2(r, t^) and 
equilibrium non-condensate mean fields, while dynamic 
shifts involve the static condensate wavefunction and 
fiuctuating non-condensate mean fields. 

In this appendix we will frequently drop space, time 
and frequency dependence in intermediate stages of the 
calculation. Time/frequency-dependent quantities can 
be interpreted directly in the time-domain; in the fre- 
quency domain, complex-conjugated terms should be 
evaluated at — w, while unstarred quantities are evalu- 
ated at +UJ. 
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1. Change in condensate energy and shape 

The change in the condensate energy and shape be- 
tween the GPE and GGPE descriptions are given by 
AA20 and A$2o(r) respectively. For consistency with 
the hnear response treatment of the fluctuations, we cal- 
culate these quantities by linearizing around the static 
solution to the ordinary GPE. Substituting Eqs. 175(1 and 
(f7S|l into Eq. lO leads to 



(GP)M$2o\ y 



-20 



Vggpe ^ 



GGPE 



(Al) 



where Vggpe denotes all terms in the static GGPE which 
are not also in the GPE [approximated with <i>2(r) ^ 
$o(i-)], i.e. 

K3GPE(r) = 2C/on(r)$o(r) + f/om^(r)$S(r) 

-/(r) + AiVof/o|$o(r)|'$o(r). (A2) 

Equation I|A1|I is solved using an expansion in the quasi- 
particle basis of Eq. 



(A3) 



where the static coefficients h^f* are given by 



^ ^_X^ L y^^^^ + < F^gpe] • ( A4) 
£j J 

The change in the condensate eigenvalue is given by 



AA20 = y"*^"" *O^GGPE 



A^ot/ol^oP ($oA$2o + *oA$ 



(A5) 



A complication which arises in the linearization lead- 
ing to Eq. l|Af |l is the fact that AA20 is not small (i.e. it 
is not proportional to the small parameter of the theory) 
except at zero temperature. The reason is that the non- 
condensate density ri(r) has a large contribution which 
is simply the ideal-gas density. In fact we can write n(r) 
in the form n(r) = nid(r) -I- nint(r), where nid(r) is the 
ideal-gas contribution and nint(r) is the part due to in- 
teractions, which is proportional to the small parameter 
of the theory. Although nid(r) is large at finite temper- 
ature, it is also roughly constant over the spatial extent 
of the condensate (this is exactly true in the Thomas- 
Fermi and homogeneous limits). As a result, the change 
in shape of the condensate A<i>2o(r) is still proportional 
to the small parameter of the theory and the linearization 
procedure is valid. 

The large size of AA20 does not cause problems else- 
where in the theory because it is cancelled straightfor- 
wardly by other terms. We can see this explicitly if 
we assume that riid(r) is in fact a constant. Its con- 
tribution to AA20 is then simply 2J7o«.id and comparison 



with Eqs. l(M|l and shows that the eS^ect of nid(r) 
in A£'4 (p, q) cancels exactly with the corresponding con- 
tribution in ^E\{j)^q). Thus only nint(r) and finite size 
effects in nid(r) are relevant for excitations and a linear 
approximation remains appropriate. 



2. Static shifts 

With the exception of the last two lines of Eq. H96|) 
for AEshapc{p,q) and the expression for AEj {p,q) of 
Eq. (|98|l . the static shifts are derived straightforwardly 
by calculating the overlap integral of the third line of 
Eq. H77|) with the quasiparticle wavefunctions for mode 
'p' using Eq. (|61|l and multiplying by — 1- This 
gives AEi{p,q), AEx{p,q), AEANo{p,q) and the first 
two lines of A£^shapo(p, <?) directly. The last two lines 
of Eq. come partly from condensate phase fluctua- 
tions and partly from the condensate norm mode, as we 
now show. 

Explicitly, the second term of Eq. H77|l is 

cf^^V^2ir,u;) = 2NoUo\<^>o\^Xg„bNo{io) 

+ [e,X,+c,Xg„]b,{uj), (A6) 

i=i±0 

where 67Vo(w) is given in Eq. I|81|) . The contribution to 
mode 'p' from the first line of Eq. (|77l) is therefore 

{huj - ep)bp{uj) - 2xpC*bNo{i^) - Xp-Ppo(w), (A7) 

where Cp and Ppo('-^) are given by Eqs. and l(7n|) 
respectively. 

The contribution from the second line of Eq. (|77|l is 

Xpjdr [P{r,cj) ~ <5A^(c.)] [^.;A<i>2o + A<i>*o] . (A8) 

In this expression, we can use the leading order result for 
dX2{uj) = dXo^Lo) given by Eq. H71|l . which removes con- 
densate phase fluctuations. Combining the result with 
Eq. l|A7p . multiplying by = 1 and regrouping we get 



Xpifi^ - ep)bp{uj) ~ Ppo{uj) = 

AP^o^)(c.)+ Y ^E^iP^l)b,{' 



(A9) 



with AP^piuj) as in Eq. lfTU2|) and AE^{p, q) deflned by 
AE^{p, q)=-c*p jdr [A$*oW, + A$20Wg] 



dr [A$20'«p + A$,nw. 



(AlO) 



AE± (p, q) depends on the change in shape of the static 
condensate wavefunction and is therefore logically in- 
cluded in the definition of Ai?shapo(?', <z), giving the last 
two lines of Eq. H9t)|l . As mentioned in the main text. 



21 



in a time-independent, perturbative calculation this con- 
tribution arises from the change in the definition of the 
subspace orthogonal to the condensate when the conden- 
sate shape changes 24]. 

The final contribution to the static shifts comes from 
Sf{r,uj) in the last line of Eq. H77|l . The corresponding 
energy shift for mode 'p' is given by 



dr [u;Sf{r,uj)+v;Sr{r,-u;)]. (All) 



AEfip) = 

Taking the expression for /(r,t) given in Eq. Hll|l and 
linearizing we obtain 

= (A12) 
where the contribution from the static non-condensate is 

Sf'^''^ = Uo y"dr'{[2|$o|''5$2 + *o'5$2]"(r,r') (A13) 

-F[2|$o|''5$* + $*2^$2]TO(r, r')}, 
while the dynamic contribution has the form 

^ l'dr'Uo\<^o\^[<i>oSn{r,r') + %67~n{r,r')]. (A14) 



All unspecified spatial coordinates are evaluated at r'. 

Using Eq. (|8()|l with h^^ — be„ — correct to this order 
gives 



^<i>2(r) - J2 "^W^?' '^'f^W = 1^ Mr)bq, (A15) 



q^±0 



while n(r, r') and fh{r,r') are obtained from the time- 
independent limit of Eqs. 1421) and H43|) . Substituting 
these results into Eqs. (|A11|) and (|A13I) leads to an ex- 
pression for the static part of the energy shift due to 
f{r,t). We are left with integrals over r' of four wave- 
functions which can be expressed in terms of the matrix 
elements Agij, Blgij and B2qij defined in Eqs. (|107() - 
(|109|l . and integrals over r of products of two functions 
which involve the factors Uij, Vij and Wij of Eqs. 
(|101|l . The final result can be written as 



AE 



(S), 



f ip) 

(S) 



(A16) 



with AEy\p,q) as in Eq. 



Expanding the field operators as in Eq. (|38|l and us- 
ing the fact that the quasiparticle operators are time- 
independent and diagonal gives 



(A18) 



where the sum is over all states and the coefficients Ri are 
defined in terms of the time-independent quasiparticle 
populations {Ni\ by 



{X^ = +1), 

0, {i = ±0). 



(A19) 



We note that the time-dependent condensate modes give 
no direct contribution in this expression because R±q = 
0. It is nonetheless convenient to include them formally 
in the summation as we have done in Eq. IjAlSp . The 
reason is that when we expand the time-dependent wave- 
functions in the time-independent basis [Eq. (|A21|l ]. the 
static condensate modes are needed to preserve orthogo- 
nality. Including them at the outset allows a symmetric 
treatment of the two summation indices which then ap- 
pear. 

We can write the time-dependent spinors in terms of 
small fluctuations around the static values given by the 
solutions of Eq. (|2H1) 



-ieit/h 



(A20) 



Expanding the fluctuations in the static basis via 



(A21) 



gives, after linearizing 



(A22) 



3. Dynamic shifts 

We now turn to the calculation of terms involving fluc- 
tuations of the non-condensate. Following Blaizot and 
Ripka we define a generalized density matrix TZ by 



nir,r',t) 



^1 (aHv'!!)) 



^ n(r, r', t) m(r, r', t) 

^m*(r,r',t) [Q* + n*] {r , r' , t) 



(A17) 



Using the time-dependent version of the orthonormal- 
ity relation of Eq. H3U|I and the relationship between pos- 
itive and negative norm spinors of Eq. H31|l we obtain the 
following useful properties of the coefficients Xij 



(A23) 



Using these results and the completeness relation of 
Eq. (|37|l . we rewrite Eq. ljA22|) in terms of a sum over 
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ij>0 



modes with positive norm only (including the state +0) 

ij>0 

X {1 + R, + Rj) X,,^, 

ij>0 

x(l + i?, + i?,)X*_^. 
J2 [A'.(r)4(r')+A'_,(r)A'!,(r'y 

X {Rj — Ri) Xji 

(A24) 




6Q*2ir,r')J ' 



where 



6Q;{r, r') = - [$o(r')(5$;(r) + n{r)S<^2{r')] . (A25) 

From this result, the formulae for 6n{r,r') and 
6m{r, r') can be read off from the first row 

6n{r,r') K(r)«;(r') +«*(r)<(r')] 

ij>0 

x{l + Ri+ Rj}X,_j 

ij>0 

X (1 + R, + Rj) X*_j 
+ [u,{r)u*{r')+v*{r)v,{r')] (A26) 

i]>0 

X {Rj — Ri) Xji, 
Srh{r,r') = ^ i [<(r)z;*(r') + z;*(r)<(r')] 

ij>Q 
ij>0 

+ Hr)v*{r')+v*{r)u.{r')] (A27) 

X {Rj — Ri) Xji. 

The dynamic terms give a contribution to the evolution 
of mode 'p' which comes from the projection of the last 
two lines of Eq. (|77jl 



dr\2UoSii[u;^o + v;^* 



(A28) 



No 



where Sf^^^ is given in Eq. ljA14p and we have used the 
fact that SNq = — Jdr Sh. Using the above equations for 
6n{r,r') and 6fh{r,r'), this becomes 



E 



{CO) 



2VM 



[1 + TV, + Nj 







(SI) 



i]>0 



Y 

ij>0 ^ 
-AE^{p,q), 



[l + N, + Nj 



N 



[Nj - N] 







(A29) 



where the coefficients xi,^^^ and fJ,^^^ are defined 

p^j^^^p'^j P''j 

in the main text in Eq. 11061) . The UV renormalization 
AE^{p,q) is given by Eq. (|105|) . This result is obtained 
straightforwardly from Eq. ljA28|) using the contribution 
to 6rn^ from the last part of Eq. (|49() . 

It remains to find an expression for the coefficients 
Xij{uj). This is achieved by solving the equation of mo- 
tion for the time-dependent Bogoliubov functions Xi{r, t) 
given in Eq. (|40l) . After linearization, this becomes 

z/i^fJA", = [£0 - e^] SX, + 6C2X,. (A30) 

with Co as in Eq. (|25|l . 6C2 is obtained by linearizing 
£2{r,t) of Eq. (Uni) using 



i/,p(r,<) = i/,p(r) + P(r,t), ^ 
$2(r,t) w $o(r) +5$2(r,t) 

g2(i) « Qo + SQ2{t), 



(A31) 



with SQ2{t) as in Eq. ljA25|) . Note that, correct to the 
order of this calculation, we have used <i>2(r) ~ ^o{^) and 
A2 ~ Ao for static condensate properties. We must retain 
the distinction between ^$o(r,t) and 5$2(r, t), however, 
for the reason given in Appendix IbI 

Substituting Eq. HA2ip into Eq. (|A30p . taking the 
Fourier transform and projecting out the required coeffi- 
cient gives 



X,j{uj) 



Xjjdr xUr)a36C2{r,u;)X,{r) 



(A32) 



applying for all states i and j. It is easily verified that 
this expression satisfies the relations of Eq. (|A23p . 

We obtain explicit expressions for Xij{uj) by writ- 
ing S£2{j^,ijj) in terms of P(r, w), (5A2(w), 5'i>2(r, t^) and 
6Q2{t)- The effect of the orthogonal projector depends 
on whether or not it acts on the condensate modes in 
the expansion. The expression for the coefficients Xij{uj) 
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therefore depends on whether or not one of the labels i 
or j corresponds to ±0. For i,j ^ ±0 and XiyXj = +1; 
we can write Xij{u}) in the form 



(A33a) 



(A33b) 



(A33c) 



In deriving these expressions the following relations be- 
tween the various coefficients we have defined are useful 



yi^) _ yy^'-)" _ y 



j,'> 

(SI)* _ ^(B2) 



(A34a) 
(A34b) 



We also note that Y^^j is symmetric with respect to i 



and J while Y^^^ is invariant with respect to permuta- 
tions of its indices. 

The form of Xij{uj) for the case that one or both of 
the indices i or j refers to ±0 is different. The expres- 
sion for Xij{uj) in this case is most easily obtained us- 
ing the fact that the time-dependent Bogoliubov func- 
tions Ui{r,t) and Vi{r,t) are orthogonal to the conden- 
sate, i.e. fdr^ l{r,t)ui{r, t) = Jdr ^2{ii^,t)vi{r,t) = 0. 
Using Eq. ijA15|l for i ^ ±0, Xi = +1 gives 

= - E ^^1^1 = (A35) 



Jdrvi 



W,,bg = X+o,-„ (A36) 



which is consistent with Eq. ljA23|) . The absence of con- 
densate phase and norm modes (at leading order) gives 



gives APp^\ll>), while AEp^\uj) comes from the part of 
Xij{cj) involving the condensate fluctuations {bq{uj)}. In 
this latter contribution, we write the part involving the 
condensate mode separately, using Eqs. HA35P and (|A36p 
which gives AEo{p,q) in Eq. (|1U4|I . We note that the 
summations in APp^\uj) of Eq. (|112|l exclude the mode 
+0 because, as in Eqs. HA35|) and (|A36|) . the terms this 
appears in can be written solely in terms of the conden- 
sate fluctuations. 

The coefficients in Eq. ljA33p are defined 

in terms of the same coefficients Yqij that appear in 
Eq. (|A29() and are defined in Eq. H1U6|I . As discussed in 
the main text, the various contributions to these coeffi- 
cients have a different physical origin in the two equa- 
tions. The derivation of Eq. IIA33p for Xij{uj) shows 
that the term involving Cg in the Yqij's conies from con- 
densate phase fluctuations [the contribution from 6X2(10) 
to SC2{r,oj)]. The corresponding term in the l^y's in 
Eq. ljA29|) comes instead from condensate number fluc- 
tuations [the term —c*Sfi/No in Eq. (|A28p ]. In addition, 
the terms involving Ci and cj in the Y's appear, in Xij{uj), 
as a consequence of fluctuations in the orthogonal projec- 
tor 6Q2- The same contributions to the Y's in Eq. HA29|I . 
come instead from Sf^^h Finally, the matrix elements 
Apij, Blpij and B2pij in the Y's describe, in AT^ (w), the 
driving of the non-condensate by the condensate fluctu- 
ation, while their appearance in Eq. I|A29|I describes the 
subsequent back action of the non-condensate on the con- 
densate. 

Finally, to calculate the condensate density fluctua- 
tions we need an expression for 6Nq(uj) — — Jdr Sn(r, lu). 
Using Eqs. HA26|I . ((TTn|) and ((TTT|l we obtain 

ij>0 
ij>0 

^ E/*X,,H[Ar,-Ar,]. (A38) 

ij>0 

This expression can be broken into two parts, correspond- 
ing to the change in n(r) due to condensate fluctuations 
and the change due to the external perturbation. The 
resulting formulae are easily obtained using the above 
results for Xij{uj). 



X 



=+n — X^c)A=+o — 



+0,i=±0 



(A37) 



Equations (|A35|I and (|A36|I show that orthogonality be- 
tween the time-dependent condensate and excited state 
wavefunctions requires bui and bvi to have a non-zero 
overlap with the static condensate modes. 

The expression for the dynamic shifts is flnally ob- 
tained by substituting Eqs. l|J33)l - l|X37|) into Eq. l|J29|l . 
The contribution to Xij (oj) from the perturbation P(r, iS) 



APPENDIX B: EFFECT OF GENERALIZED GPE 
WAVEFUNCTION IN NON-CONDENSATE 
EQUATIONS 

A subtle issue in the theoretical formalism concerns the 
appearance of the generalized GPE wavefunction <I>2(r,t) 
rather than <I>o(r,i) in the matrix £2(1", i) which defines 
the evolution of the non-condensate [c.f. Eqs. H19|) - (|2I|) ]. 
In a self-consistent treatment of the non-condensate dy- 
namics £2(1", would also contain terms involving ri(r, i) 
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and m(r, t) describing the action of non-condensate fluc- 
tuations on themselves. Such a treatment has been given 
in 1^ using a symmetry-breaking approach, and also in- 
cludes additional kinetic terms in both the condensate 
and non-condensate equations of motion which are re- 
quired for consistency. In this paper we treat the coupling 
between the condensate and non-condensate perturba- 
tively so that the definition of £2(r, t) given in Eqs. H19|l- 
(I21|l is sufficient for our purposes. 



However, the use of $2(1", rather than $o(r, in 
the equation of motion for the non-condensate produces 
some theoretical difficulties. For example, orthogonality 
between the condensate and excited states in the time- 
dependent case is only preserved to the order of the the- 
ory and is not exact. We stress that this is not a problem 
for the static basis we use to expand wavefunctions be- 
cause this is defined by Eq. (^5)1 where the replacement 
'I'2(i") — > <fo(r) is explicitly used, but the issue does arise 
when we solve for the dynamic fluctuations of the non- 
condensate, as in Appendix I A 31 Nonetheless, it is essen- 
tial to use the generalized wavefunction in the dynamic 
case if we wish to obtain sensible expressions for energy 
shifts, as we demonstrate below. This should be valid 
within the linearized, perturbative approach of interest 
here, but a self-consistent treatment will be needed for 
higher order calculations or for situations where the con- 
densate is perturbed beyond the linear regime. 



For the purposes of this appendix, we distinguish be- 
tween the condensate excitation coefficients {biltu)} ob- 
tained from the GPE and GGPE theories by denoting 
the former by {bf\Lu)} and the later by {bf\uj)}. The 
evolution of the non-condensate introduces the dynamic 
energy shift AeI^\uj) and the change in the excitation 
matrix element APpq\lu) into the theory. Of these, 

only AEp^\uj) depends on the treatment of the conden- 
sate dynamics. This quantity is written in Eq. (|84|l and 
depends on the expansion coefficients, which should be 
taken to be {bf\uj)} rather than {b^l^\uj)} because of the 
use of $2(1", i) in the non-condensate evolution. Special- 
izing to the case that only one quasiparticlc mode is ex- 
cited, this means that when we solve Eq. (|82|l AEp^\lo) 
must be taken over to the left-hand-side of the equa- 
tion (which contains bp ) and hence appears in the de- 
nominator of the resulting expression for Qp{ijj) given in 
Eq. (|S5|) . If instead we were to use $o(r,t) to evolve 
the non-condensate then AEp^'^ (w) would be defined in 

terms of the previously calculated coefficients {bf'^u;)} 
given in Eq. (|67|l . For consistency, we would also have to 
treat the static shifts of Eq. (|83ll in the same approxima- 
tion (to avoid infra-red divergences). The result is that 
the self-energy Ep(u;) of Eq. H89|) would appear in the nu- 
merator of Gp{uj) rather than the denominator, and we 



would obtain 



6(2) 
P 



huj — e„ 



Sp(^) 

huJ — €„ 



(Bl) 



Comparing this with Eqs. 
sponds to writing 



l|HSll-l|HHI)i we see that it corre- 



1 



(B2) 



{huj -~ ep) [1 - Y.p{uj) / {huj - ep)] ' 

and expanding the denominator to first order in 
Y,p{uj) / [huj — Ep). This is certainly not a valid approx- 
imation in the interesting region ftw ~ e^, and does not 
lead to a shift in the peak of Qp{uj) by an amount cor- 
responding to Y,p{Lo). It is therefore necessary to include 
the condensate dynamics self-consistently in the theory 
to obtain sensible expressions for energy shifts and this is 
the reason why $2 (r, t) appears in the equation of motion 
for the non-condensate. 



APPENDIX C: INTRODUCTION OF AN 
IMAGINARY FREQUENCY IN THE 
RESOLVENT 

The expressions for the dynamic contributions to the 
energy shift AEp^'^ (w) and excitation matrix element 

APpq \uj) contain energy denominators which may van- 
ish or become very small if a particular collision pro- 
cess is energetically allowed. To obtain finite quanti- 
ties, essential for numerical work, it is therefore neces- 
sary to include a small imaginary part in the frequency 
via uj ^ uj + i^ [33, ■ For a homogeneous gas, where an 
excitation couples to a continuum of decay channels, we 
can take the limit 7^0"*", but for trapped condensates 
the discrete nature of the states involved means that we 
must work with a finite (but small) value of 7. It is there- 
fore necessary to justify its inclusion more carefully and 
to predict its approximate numerical value. 

We can justify the inclusion of 7 by considering the 
finite experimental observation time. In fact, the coeffi- 
cients {bi{uj)} with 7 = do not correspond to any ex- 
periments because they assume that the GPE describes 
the measured evolution of the condensate for all time. 
This of course is not the case and in any real experi- 
ment the condensate must first be created and allowed 
to settle into its equilibrium state. It is then excited 
by modulation of the trapping potentials and the re- 
sulting disturbance is observed for some finite length of 
time. Thus the quantity we are actually interested in 
experimentally is not (5$(r, i), but rather something like 
0(O)5$(r, t)9(robs — t), where Tobs is the experimental 
observation time and Q{t) is the unit step function. Mod- 
elling 0(Tobs — by the decaying exponential e~*/'^°^' we 
therefore want to calculate 



5$(r,i) = (5$(r,t)e~'''e(0), 



(CI) 
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with 7 = 1/Tobs and 5$(r, 0) = 0. The Fourier transform 
of this quantity amounts to putting a; ^ w + 17 in the 
coefficients {bi{u))}. This replacement can be neglected 
in the perturbation P{r, co) because this is not divergent 
and will generally have an intrinsic width greater than 7 
in any case. 

If more accurate modelling of the experimental reso- 
lution is required, numerical results for the condensate 
response in the frequency domain can be transformed to 
the time-domain at the end of the calculation. The factor 
e"'''' can then be removed and some other weight function 
introduced as appropriate. Of course, when calculating 
decay rates by fitting Lorentzians to the calculated spec- 
tra, the explicit contribution from 7 in the denominator 
of Gp{uj + ij) in Eq. (|S5|) must be subtracted. The re- 
maining dependence of the decay rate and energy shift on 
7 comes from the self-energy Ep(aj + 17) of Eq. (|89|l . In 
practice this dependence is usually weak and numerical 
results do not depend significantly on 7 for experimen- 
tally relevant values (typically of order a few times 10"^ 
of the trapping frequency). 

Of course, there may be sources of broadening other 
than experimental resolution which can conveniently be 
included in 7, such as collisions with background gas 
atoms for example. In principle, we could also model non- 
condensate collisions as an additional source of broaden- 
ing of the quasiparticle levels. This would only affect 
the value of 7 in the dynamic terms AEp^"^ (lo -f 17) and 

APp,^' (to -\- ij) and would involve having two values of 7 
in the theory, an explicit contribution 71 in the denomi- 
nator of Gp and a larger value 72 in the dynamic terms. 
If this had a significant effect, however, such processes 
would have to be described more rigorously via a higher 
order extension of the theory. 



{Nq) — {No)"^, where the averages are calculated in 
the quasiparticle basis. For a system with a fixed to- 
tal number of atoms wc have Var(A^o) = Var(A^nc) = 
(Nl^) - (iV„e)2, where N^c = /dr At(r, i)A(r, i) is the 
non-condensate number operator. Dropping the time la- 
bel, (N^c) is therefore defined by 



//.'rdV (A. WAWA.„')A(r')). (Dl) 



This can be calculated in equilibrium using Wick's the- 
orem which gives 



At(r)A(r)A'f (r')A(r')) = n(r)n(r') + |m(r, r')P (D2) 
+ n(r',r) [n(r, r') + Q(r, r')] . 

We therefore have 
Var(Ar„^) = Nnc + ffd^rd^r' |n(r,r')P + |m(r,r')P, 



(D3) 

where TVnc = (-^nc) = /drn(r). After some algebra this 
reduces to 



Var(7V„c) = 



ij>0 



ij>0 



(D4) 



APPENDIX D: VARIANCE IN THE 
CONDENSATE POPULATION 

To calculate AA^o of Eq. H12(l . we need an expression 
for the variance in the condensate population Var(A^o) = 



where lij, Jij, are defined in Eqs. H11Q|1 and Hlll(l while 
Uij and Vij are defined in Eqs. (|^ and H1UU|) . In the 
homogeneous limit this expression agrees with the result 
obtained in 18(1. 
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